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Three  Dimensional  Transient  Analysis  of  Microstrip  Circmts  in  Miiltilavered  Anisotropic  Media 


Under  the  sponsorship  of  the  ONR  contract  N00014-90-J-1002  we  have  published 
(18)  refereed  journal  and  conference  papers. 

In  high-speed  digital  circuits,  high  frequency  phenomena  affects  the  characteristics 
of  the  interconnections.  Physical  discontinuities  or  nonuniformities  in  the  connections  may 
cause  severe  reflections  when  they  can  no  longer  be  considered  as  conducting  wires,  but 
behave  as  transmission  lines  and/or  waveguides.  In  multilayered  digital  circuits,  vias  con¬ 
stitute  one  of  the  most  commonly-used  class  of  interconnects.  Vias  are  not  good  czirriers 
of  high-speed  signals.  They  cause  signail  distortion  and  reflections  as  well  as  severe  degra¬ 
dation  in  the  high  frequency  components.  The  analysis  smd  modeling  of  a  single  via  and 
some  quasi-static  or  quasi-TEM  analyses  of  single  via  configurations  have  been  carried  out 
previously.  To  date,  the  analysis  of  coupled  noise  between  adjacent  vias  hsis  received  very 
little  attention.  The  major  re^lsons  include  the  complexity  of  multi- via  structures  and  the 
difficvilty  of  modeling  and  analyzing  them  accurately  over  a  broad  frequency  range.  The 
electromagnetic  coupling  between  two  adjacent  vias  in  a  multilayered  integrated  circuit  is 
analyzed  by  means  of  equivalent  magnetic  frill  array  models  incorporated  with  the  even- 
and  odd-mode  approach.  Closed-form  expressions  for  the  coupled  noise  on  the  passive  via 
are  derived.  The  coupling  responses  in  the  frequency  domain  and  crosstalk  waveforms  in 
the  time  domain  for  some  multilayered  via  structures  are  calculated  based  on  these  for¬ 
mulas.  A  4-layer  experimental  model  is  constructed  and  measurements  are  taken  for  the 
trzmsmission,  reflection,  and  coupling  responses.  The  me2isurements  show  good  agreement 
with  the  calculated  results  over  a  frequency  range  of  up  to  18  GHz. 

High  frequency  digital  signals  associated  with  computers  and  communications  equip- 
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ment  have  the  potential  for  electromagnetic  interference  (EMI).  The  need  for  comphance 
with  the  radiated  emissions  limits  established  by  FCC  and  other  regulatory  agencies  for 
such  equipment  have  led  to  increased  interest  in  understanding  and  minimizing  EM  radia¬ 
tion.  As  clock  speeds  increase  beyond  100  MHz  and  edge-rates  fall  into  the  sub-nanosecond 
range,  the  EMC  engineer  is  faced  with  FCC  tests  for  frequencies  beyond  1  GHz  and  the 
prospect  of  increased  energy  for  frequencies  in  the  hundreds  of  megahertz  range.  The 
modules-on-backplane  configuration  is  a  common  assembly  scheme  for  computers.  Em¬ 
bedded  in  these  modules  and  backplane  are  reference  planes  which  could  be  as  large  as 
tens  of  centimeters.  The  reference  planes  configuration  is  of  interest  because  the  resonant 
frequencies  begin  in  the  low  hundreds  of  megahertz,  given  typical  dimensions.  The  voltage 
levels  on  such  reference  planes  are  subject  to  noise  fiuctuations  (reference  noise)  originat¬ 
ing  from  disturbances  such  as  inductive  discontinuities  at  connector  pins.  The  increase  in 
power  requirements  and  in  the  number  of  simultaneous  switching  circuits  is  likely  to  result 
in  stronger  noise  sources.  A  simplified  model  is  used  to  analyze  the  modules-on-backplane 
configuration.  Conducting  planes  model  the  modules  and  the  backplane.  A  generator, 
which  is  positioned  at  the  module  and  backplane  interface,  models  the  inductive  voltage, 
as  a  noise  source,  developed  across  the  module-backplane  connector.  The  finite- difference 
time-domain  (FD-TD)  technique,  which  is  based  on  the  discretization  of  the  Maxwell’s 
equations,  is  employed  in  this  analysis.  Nodal  electric  and  magnetic  fields  are  computed 
on  a  rectangiilar  grid  within  a  computationsd  domain  subject  to  appropriate  boundary  and 
initial  conditions.  An  absorbing  boundary  condition  is  enforced  on  the  outer  boundary  of 
the  computational  domain  and  approximately  simulates  unbounded  space.  In  the  imple¬ 
mentation  of  the  FD-TD  technique,  the  areas  of  concern  addressed  include  the  modeling 
of  the  excitation  as  a  coaxial  line  probe  and  the  use  of  Prony’s  method  to  obtain  late 
time  responses  of  the  fields  for  highly  resonant  configurations.  Radiation  impedances  are 
obtained  from  the  voltages  and  currents  on  the  probe  and  radiation  patterns  are  calculated 
using  Huygens’  principle.  Experimental  verification  is  made  on  physicsd  structures  com- 
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patible  with  simulation  models.  The  focus  is  on  correlating  computed  resonant  frequencies 
and  radiation  impedances  with  those  deduced  from  ^-parameter  measurements  taken  on 
a  network  analyzer.  Radiation  impedances  as  functions  of  frequency  are  presented  for 
various  configurations.  Such  variations  include  changing  the  number  of  modiiles  and  the 
dimensions  and  spacings  of  modules  and  backplane  and  the  addition  of  conducting  panels 
to  model  the  shielding  or  enclosure  environment.  This  study  provides  basic  information 
pert^lining  to  the  acceptability  of  noise  levels  in  the  interconnect  environment  for  EMI- 
related  considerations.  Through  the  examination  of  a  variety  of  configurations  it  will  also 
affect  system-level  packaging  and  assembly  decisions. 

Millimeter  and  subnullimeter  transmission  structures  are  dominated  by  dielectric 
wavegmdes.  At  these  frequency  bands,  conventional  microstrip  trai^smission  structures 
such  as  microstrip  line,  coplanar  line,  slot  line  and  finline  all  suffer  from  severe  ohmic  and 
radiation  losses.  Analytical  techniques  for  studying  dielectric  waveguide  structures  include 
the  full-wave  integral  equation  formulations  (both  spectral  and  spatial  domain)  and  the  fi¬ 
nite  difference  time-domain  (FDTD)  method.  Compared  to  the  integral  equation  methods, 
the  FDTD  method  is  more  straightforward  to  formulate  and  provides  an  alternative  path 
to  analyze  both  time  and  frequency  domain  response  in  the  presence  of  complex  structures 
such  as  discontinuities.  In  this  study,  electromagnetic  pulse  propagation  in  single  and 
coupled  dielectric  rib  wavegvdde  embedded  in  a  three-layer  stratified  dielectric  medium  is 
simulated  and  analyzed  using  the  FDTD  technique.  Time  domain  transient  response  of  a 
modulated  gaussian  pulse  is  observed  in  the  simulation.  We  first  examine  pulse  dispersion, 
attenuation  and  coupling  characteristics  for  the  transmission  line  structure.  Results  are 
found  to  be  in  good  agreement  with  those  obtained  from  integral  equation  methods.  We 
also  extend  the  study  to  an  optical  rib  waveguide  with  a  bend  discontinuity.  The  FDTD 
simulation  is  used  to  predict  the  electromagnetic  leakage  attributed  to  the  discontinuity. 
One  of  the  problems  facing  the  FDTD  technique  for  the  dielectric  waveguide  problem  is 
the  proper  absorbing  boundary  condition.  Existing  boundary  conditions  such  eis  the  Mur, 
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Enquist-Majda  condition  work  best  when  wave  incident  angle  is  normal  to  the  boundary; 
which  is  generally  not  the  case  for  dielectric  waveguides.  A  comparison  between  different 
conditions  will  also  be  presented. 

The  use  of  equipment  enclosures  to  shield  against  emissions  and  external  energy  is 
a  common  practice.  Efforts  range  from  metal  impregnated  plastics  to  thick  steel  panels 
incorporating  waveguide  design  in  air  vents.  The  choice  of  complexity  depends  on  system 
performance  and  application.  The  question  has  often  been  raised  as  to  whether  high  field 
intensities  within  such  enclosures  will  affect  eqvupment  performance.  Clearly  this  is  in¬ 
creasingly  pertinent  in  cases  of  heavy  shielding  where  resonator-like  structures  with  high 
quality  factors  (Q)  result.  Moreover,  in  these  cases,  any  energy  leakage  may  be  highly 
frequency  selective  with  substantial  associated  field  strength.  The  primary  motivation  for 
studying  such  structures  is  that  computers  are  now  operating  at  faster  speeds  and  consum¬ 
ing  more  power  resulting  in  significantly  increased  levels  of  power  at  higher  frequencies. 
Also  typical  metallic  enclosure  configurations  have  dimensions  which  are  electrically  reso¬ 
nant  at  frequencies  in  the  hundreds  of  megahertz  range.  We  focus  on  the  use  of  resistive 
material  within  the  enclosure  walls  for  absorbing  some  of  the  contained  RF  energy,  thereby 
avoiding  high  Q  enclosures  and  reducing  the  likelihood  of  significant  field  buildup  within 
the  enclosure.  The  energy  leakage  from  perforations  on  the  walls  is  also  being  considered. 
The  radiation  properties  of  these  structures  are  observed  by  varying  the  number,  size  and 
position  of  the  holes  in  the  enclosure.  The  finite-difference  time-domain  (FD-TD)  tech¬ 
nique  is  used  in  analyzing  this  problem  and  the  metrics  are  total  radiated  power  and  field 
patterns  within  the  enclosure.  The  excitation  sources  are  line  current  source  for  the  two- 
dimensional  models  and  dipoles  for  the  three-dimensional  cases.  A  Gaussian  amplitude 
is  assumed  in  order  to  obtain  multifrequency  characteristics  and  to  exhibit  the  resonance 
properties  of  the  model.  In  addition  to  the  radiated  power  through  the  apertures,  the 
power  absorbed  in  the  resistive  material  is  also  calculated. 
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A  new  method  based  on  the  wave  transnussion  matrix  is  presented  to  deal  with 
nonuniform  multiconductor  line  system.  This  method  has  not  only  the  ability  of  han¬ 
dling  nonuniformity,  but  also  the  advantage  of  handling  nonlinearity  efficiently.  The  wave 
transmission  matrix  is  employed  to  characterize  the  multiconductor  transmission  lines  and 
the  telegrapher  equations,  therefore,  can  easily  be  turned  into  algebraic  equations.  The 
nonuniform  transmission  line  is  treated  iis  many  small  segments  of  uniform  transmission 
lines  connected  in  cascade.  The  formulation  for  this  method  is  described  in  detail.  Calcu¬ 
lation  approaches  are  discussed.  Numerical  results  are  presented  and  comparison  of  them 
with  previous  work  confirms  the  validity  of  the  method. 

A  macroscopic  model  is  proposed  for  nonlinear  electromagnetic  phenomena  in  super¬ 
conductors.  Nonlinear  constitutive  relations  are  derived  by  modifying  the  linear  London’s 
equations.  The  superelectron  number  density  as  a  function  of  applied  macroscopic  current 
density,  ns(J)^  is  derived  from  a  distribution  of  electron  velocities  at  a  certain  temperature 
T.  At  temperature  T  ^  OK,  the  function  ns{J)  has  a  smooth  variation  near  the  macro¬ 
scopic  critical  current  density  Jc.  Agreement  has  been  found  between  this  ns{J,T)  model 
and  the  temperature  dependence  of  in  the  two-fluid  model.  The  nonlinear  conductivities 
o‘s{J)  and  <rn{J)  axe  obtained  from  the  London’s  equation  with  the  modified  ns{J)  func¬ 
tion.  Nonlinear  resistance  R{I),  kinetic  inductance  L]g{I)  and  surface  impedance  Zs{I)  in 
thin  wire,  slab,  and  strip  geometries  are  calculated. 

The  finite  difference-time  dommn  (FD-TD)  technique  is  applied  to  the  solution  of 
Maxwell’s  equations.  A  computer  program,  which  can  be  used  to  simiilate  and  study 
numerous  electromagnetic  phenomena,  is  developed  and  implemented  on  a  386-class  per¬ 
sonal  computer.  The  FD-TD  technique  is  a  useful  tool  for  students  in  electromagnetics. 
The  technique  is  flexible  and  can  be  applied  to  many  basic  EM  scattering  and  radiation 
problems.  Because  field  solutions  are  found  as  a  function  of  time,  visualization  of  the 
propagation  of  the  EM  fields  is  possible.  The  FD-TD  technique  is  implemented  for  a 
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two-dimensional  rectangular  grid  in  conjunction  with  a  second-order  absorbing  boundary 
condition.  Both  £1-  and  H-field  polarizations  are  analyzed.  Finite  objects  consisting  of 
dielectric,  magnetic  and  conducting  materials,  and  perfectly  conducting  infinite  ground 
planes  sire  modeled.  Plsme  wave  and  line  current  sources  are  implemented.  In  addition 
to  the  capability  of  animating  the  propagation  of  the  EM  fields,  radiation  and  scattering 
patterns  can  be  generated. 

A  methodology  developed  to  handle  dispersive  materials  in  the  time  domain  is 
extended  to  model  the  dispersive  characteristics  of  the  impedance  boundary  condition 
used  for  a  thin  layer  coating  over  perfect  conductors.  The  impedance  boundary  condition 
is  first  approximated  as  a  rational  function  of  frequency.  This  rational  function  is  then 
transformed  to  a  time  domain  equation,  resulting  in  a  partial  differential  equation  in 
space  and  time.  Discretization  of  the  time  domain  model  to  efficiently  handle  the  thin 
layer  coating  is  presented  in  the  context  of  the  finite-difference  time-domain  (FD-TD) 
technique.  The  methodology  is  verified  by  solving  a  one-dimensional  problem  using  the 
FD-TD  technique  and  comparing  with  the  analytical  results. 

To  understand  the  physical  meaning  of  rational  reflection  coefficients  in  inverse¬ 
scattering  theory  for  optical  waveguide  design,  we  studied  the  relationship  between  the 
poles  of  the  transverse  reflection  coefficient  and  the  modes  in  inhomogeneous  dielectrics. 
By  using  a  stratified-medium  formulation  we  showed  that  these  poles  of  the  spectral  re¬ 
flection  coefficient  satisfy  the  same  equation  as  the  guid^mce  condition  in  inhomogeneous 
waveguides.  Therefore,  in  terms  of  wave  numbers,  the  poles  are  the  S2mie  as  the  discrete 
modes  in  the  waveguide.  The  radiation  modes  have  continuous  real  values  of  tremsverse 
wave  numbers  and  are  represented  by  the  branch  cut  on  the  complex  plane.  Based  on 
these  results,  applications  of  the  Gel’fand-Levitan-Marchenko  theory  to  opticzd  waveguide 
synthesis  with  the  rational  function  representation  of  the  transverse  reflection  coefficient 
are  discussed. 
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The  coupled'Wave  theory  is  generalized  to  analyze  the  diffraction  of  waves  by  chiral 
gratings  for  arbitrary  soxgles  of  incidence  and  poleirizations.  Numerical  results  for  the 
Stokes  parameters  of  diffracted  Floquet  modes  versus  the  thickness  of  chiral  gratings  with 
various  chiralities  are  calculated.  Both  horizontal  and  vertical  incidences  are  considered  for 
illustration.  The  diffracted  waves  from  chiral  gratings  are  in  general  elliptically  polarized; 
and  in  some  particular  instances,  it  is  possible  for  chiral  gratings  to  convert  a  linearly 
polarized  incident  field  into  two  nearly  circularly  polarized  Floquet  modes  propagating  in 
different  directions. 

A  general  spectral  domain  formulation  to  the  problem  of  radiation  of  arbitrary 
distribution  of  sources  embedded  in  a  horizontally  stratified  arbitrary  magnetized  linear 
plasma  is  developed.  The  fields  are  obtained  in  terms  of  electric  and  magnetic  type  dyadic 
Green’s  functions.  The  formulation  is  considerably  simplified  by  using  the  kDB  system 
of  coordinates  in  conjunction  with  the  Fourier  transform.  The  distributional  singular 
behavior  of  the  various  dyadic  Green’s  functions  in  the  source  region  is  investigated  and 
taken  into  account  by  extracting  the  delta  function  singularities.  Finally,  the  fields  in  any 
arbitrary  layer  are  obtained  in  terms  of  appropriately  defined  global  upwzird  and  downward 
reflection  and  transmission  matrices. 

We  investigated  a  method  for  the  calculation  of  the  current  distribution,  resistsmce, 
and  inductance  matrices  for  a  system  of  coupled  superconducting  tr2msmission  lines  having 
finite  rectangular  cross  section.  These  calc\ilations  allow  accurate  characterization  of  both 
high-Tc  and  low-Tc  superconducting  strip  transmission  lines.  For  a  single  stripline  geome¬ 
try  with  finite  ground  planes,  the  current  distribution,  resistance,  inductance,  and  kinetic 
inductance  are  calculated  as  a  function  of  the  penetration  depth  for  various  film  thickness. 
These  calculations  are  then  used  to  determine  the  penetration  depth  for  Nb,  NbN,  and 
Y Ba2Cu^O'j_x  superconducting  thin  films  from  the  measured  temperature  dependence 
of  the  resonant  frequency  of  a  stripline  resonator.  The  calculations  are  also  used  to  convert 
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measured  temperature  dependence  of  the  quality  factor  to  the  intrinsic  surface  resistance 
as  a  function  of  temperature  for  a  Nb  stripline  resonator. 

The  electromagnetic  radiation  from  a  VLSI  chip  package  and  heatsink  structure 
is  analysed  by  means  of  the  finite-difference  time-domain  (FD-TD)  method.  The  FD- 
TD  algorithm  implemented  incorporates  a  multi-zone  gridding  scheme  to  accommodate 
fine  grid  cells  in  the  vicinity  of  the  heatsink  and  package  cavity  and  sparse  gridding  in 
the  remainder  of  the  computational  domain.  The  issues  pertaining  to  the  effects  of  the 
heatsink  in  influencing  the  overall  radiating  capacity  of  the  configuration  zure  addressed. 
Analyses  are  facilitated  by  using  simplified  heatsink  models  and  by  using  dipole  elements  as 
sources  of  electromagnetic  energy  to  model  the  VLSI  chip.  The  potential  for  enhancement 
of  spurious  emissions  by  the  heatsink  structure  is  examined.  For  heatsinks  of  typical 
dimensions,  resonance  is  possible  within  the  low  gigahertz  frequency  range. 

Because  the  effects  of  diffraction  during  proximity-print  x-ray  lithography  are  of 
critical  importance,  a  number  of  previous  researchers  have  attempted  to  calculate  the 
diffraction  patterns  and  minimum  achievable  feature  sizes  as  a  function  of  wavelength  and 
gap.  Work  to  date  has  assumed  that  scalar  diffraction  theory  is  applicable-as  calculated,  for 
example,  by  the  Rayleigh-Sommerfeldformulation-and  that  Kirchhoff  boundary  conditions 
can  be  applied.  Kirchhoff  boundary  conditions  assume  that  the  fields  (amplitude  and 
phase)  are  constant  in  the  open  regions  between  absorbers,  and  a  different  constant  in 
regions  just  under  the  absorbers  (i.e.,  that  there  are  no  bringing  fields).  An  x-ray  absorber 
is,  however,  best  described  as  a  lossy  dielectric  that  is  tens  or  hundreds  of  wavelengths 
tall,  and  hence  Kirchhoff  boundary  conditions  are  unsuitable.  We  have  instead  used  two 
numerical  techniques  to  calculate  accurate  diffracted  fields  from  gold  absorbers  for  two 
cases:  a  30  nm-wide  line  at  A  =  4.5  nm,  and  a  100  nm-wide  line  at  A  =  1.3  nm.  We 
show  that  the  use  of  Kirchhoff  boundary  conditions  introduces  unphysically  high  spatial 
frequencies  into  the  diffracted  fields.  The  suppression  of  these  frequencies-which  occurs 
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naturally  without  the  need  to  introduce  an  extended  source  or  broad  spectrum-improves 
exposure  latitude  for  mask  features  near  0.1  fim  and  below. 

In  order  to  understand  the  physical  meaning  of  rational  reflection  coefficients  in 
one-dimensioned  inverse  scattering  theory  for  optical  waveguide  design,  we  have  studied 
the  relation  between  the  poles  of  the  transverse  reflection  coefficient  and  the  modes  in 
inhomogeneous  dielectrics.  By  using  a  stratified  medium  model  it  is  shown  that  these 
poles  of  the  reflection  coefficient  have  a  one-to-one  correspondence  to  the  discrete  modes, 
which  are  the  guided  and  leaky  modes.  The  radiation  modes  have  continuous  real  values  of 
transverse  wave  numbers  and  are  not  represented  by  the  poles  of  the  reflection  coefficient. 
Based  on  these  results,  applications  of  the  Gel’fand-Levitan-Marchenko  theory  to  optical 
wavegmde  synthesis  with  the  rational  function  representation  of  the  transverse  reflection 
coefficient  are  investigated. 

In  compact  modules  of  high  performance  computers,  signal  transmission  lines  be¬ 
tween  integrated  circuit  chips  are  embedded  in  mtiltilayered  dielectric  medium.  These 
signal  lines  are  usually  placed  in  different  layers  and  run  perpendiculu  to  each  other.  The 
interaction  between  the  orthogonal  crossing  lines  and  the  signal  line  affects  its  propagation 
characteristics  and  may  cause  considerable  signal  distortion. 

The  interaction  of  a  pair  of  crossing  lines  in  isotropic  medium  has  been  studied  using 
a  time-domain  approach,  where  coupling  is  described  qualitatively.  This  method  becomes 
computationally  expensive  when  the  number  of  crossing  lines  increases.  With  many  identi¬ 
cal  crossing  strips  uniformly  distributed  above  the  signal  line,  the  transmission  properties 
are  characterized  by  stopbands  due  to  the  periodicity  of  the  structure.  Periodic  struc¬ 
ture  have  been  investigated  using  frequency-domain  methods.  Periodically  nonuniform 
microstrip  lines  in  an  enclosure  are  analyzed  on  the  basis  of  a  numerical  field  calculation. 
A  technique  based  on  the  network-analytical  formulism  of  electromagnetic  fields  has  been 
used  to  analyze  striplines  and  finlines  with  periodic  stubs.  The  propagation  characteristics 
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of  waves  along  a  periodic  array  of  parcdlel  signal  lines  in  a  multilayered  isotropic  struc¬ 
ture  in  the  presence  of  a  periodically  perforated  ground  pl2me  and  that  in  a  mesh-plane 
environment  have  been  studied.  More  recently,  the  effect  of  the  geometrical  properties  on 
the  propagation  characteristics  of  strip  lines  with  periodic  crossing  strips  embedded  in  a 
shielded  one-layer  isotropic  medium  have  been  investigated.  Both  open  and  closed  mul¬ 
tilayered  uniaxially  anisotropic  structures  are  considered.  A  full-wave  analysis  is  used  to 
study  the  propagation  characteristics  of  a  microstrip  line  in  the  presence  of  crossing  strips. 
The  signal  line  and  the  crossing  strips  are  assumed  to  be  located  in  two  arbitrary  layers  of 
a  stratified  uniaxially  anisotropic  medium.  An  integral  equation  formulation  using  dyadic 
Green’s  functions  in  the  periodically  loaded  structure  is  derived.  Galerkin’s  method  is 
then  used  to  obtain  the  eigenvalue  equation  for  the  propagation  constant.  The  effects  of 
anisotropy  on  the  stopband  properties  are  investigated.  Numericzd  results  for  open  and 
shielded  three-layer  unizixially  anisotropic  media  are  presented. 

For  microwave  integrated  circuit  applications,  the  characteristics  of  interconnects 
have  been  investigated  for  the  propagation  modes,  time  response,  crosstalk,  coupling, 
delay,  etc.  In  these  analyses,  it  is  assumed  that  quasi-TEM  modes  are  guided  along  the 
multiconductor  transmission  lines.  The  analysis  were  performed  for  arbitrary  number  of 
transmission  lines  where  the  load  and  the  source  conditions  were  presented  in  terms  of  the 
modal  reflection  and  transmission  coefficient  matrices. 

To  perform  the  quasi-TEM  analysis,  the  capacitatnce  matrix  for  the  multiconductor 
transmission  line  hais  to  be  obtained  first.  Both  the  spectral  and  the  spatial  domain 
methods  have  been  proposed  to  calculate  the  capacitance  matrix.  In  the  spectral  domain 
methods,  two  side  walls  are  used  to  enclose  the  whole  transmission  line  structure,  and  the 
thickness  of  the  strip  lines  has  not  been  considered.  In  using  the  spatial  domain  method, 
the  structure  has  to  be  truncated  to  a  finite  extent  to  make  the  numerical  implementation 
feasible.  The  infinite  extent  of  the  structure  was  also  incorporated,  but  only  a  two-layer 
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medium  was  considered. 


In  practical  microwave  integrated  circuits,  the  dielectric  loss  due  to  the  substrate 
and  the  conductor  loss  due  to  the  metallic  strips  are  zJso  studied  in  the  analysis  of  circuit 
performances. 

Based  on  the  scalar  Green’s  function,  a  set  of  coupled  integral  equations  is  obtained 
for  the  charge  distribution  on  the  strip  surfaces.  Pulse  basis  functions  and  a  point-matching 
scheme  is  used  to  solve  numerically  the  set  of  integral  equations  for  the  ch2n:ge  distribution, 
and  hence  the  capacitance  matrix.  The  duality  between  the  electrostatic  formiilation  and 
the  magnetostatic  one  is  applied  to  calculate  the  inductance  matrix.  The  conductance 
matrix  is  obtained  by  using  the  duality  between  the  electrostatic  problem  and  the  current 
field  problem.  A  perturbation  method  is  used  to  calculate  the  resistance  matrix. 

Finally,  a  transmission  line  analysis  is  derived  to  obtain  the  transfer  matrix  for  multi- 
conductor  uniform  lines,  which  significantly  reduces  the  effort  in  treating  the  load  and  the 
source  conditions.  Transient  responses  are  obtained  by  using  the  Fourier  transform.  The 
results  for  two  coupled  lines  are  obtained. 

With  the  ever  increasing  speed  and  density  of  modern  integrated  circuits,  the  need 
for  electromagnetic  wave  analysis  of  phenomena  such  as  the  propagation  of  transient  sig¬ 
nals,  especially  the  distortion  of  signal  piilses,  becomes  crucial.  One  of  the  most  important 
causes  of  pulse  distortion  is  the  frequency  dependence  of  conductor  loss,  which  is  caused 
by  the  “skin  effect”,  and  which  can  be  incorporated  into  the  circmt  models  for  transmis¬ 
sion  lines  as  firequency-dependent  resistance  and  inductance  per  unit  length.  Efficient  and 
accurate  algorithms  for  C2ilculating  these  parameters  are  increasingly  importeint. 

We  have  developed  a  hybrid  cross-section  finite  element /coupled  integral  equation 
method.  The  technique  is  a  combination  of  a  cross-section  finite  element  method,  which  is 
best  for  high  frequencies.  An  interpolation  between  the  results  of  these  two  methods  gives 
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very  good  results  over  the  entire  frequency  range,  even  when  few  basis  functions  are  used. 


In  the  cross-section  method,  we  divide  each  conductor  into  triangulsu'  patches  and 
choose  one  of  the  patches  from  the  return  conductor  to  be  our  reference.  We  then  calculate 
the  resistance  and  inductance  matrices  for  the  patches.  Using  two  conditions  on  the  system, 
that  the  total  current  in  each  vrire  is  the  sum  of  the  currents  in  the  patches,  and  that  the 
voltage  on  each  patch  in  a  wire  must  be  the  same  (no  transverse  currents),  we  can  reduce 
the  matrices  for  the  patches  to  the  matrices  for  the  wires.  In  the  Weeks  method,  the 
patches  are  rectangles,  and  the  quadruple  integral  is  done  quite  easily  in  closed  form. 
However,  it  is  also  possible  to  evaluate  the  quadruple  integral  in  closed  form  for  triangular 
patches,  although  the  mathematics  leading  to  this  result  is  quite  involved,  and  the  final 
form  of  the  answer  is  complicated.  We  therefore  use  triangular  patches  as  the  most  flexible 
means  of  modelling  conductors  with  arbitrary  cross-sections;  polygons  are  covered  exactly, 
and  we  are  able  to  model  quite  closely  other  shapes,  such  as  circles. 

As  frequency  increases,  the  need  to  keep  the  uniform  current  approximation  valid  in 
the  patches  requires  either  the  addition  of  many  more  patches  as  the  skin  depth  decreases, 
or  a  redistribution  of  the  existing  patches  to  the  surface,  where  the  current  is.  However, 
changing  the  distribution  of  patches  makes  it  necessary  to  recalculate  the  resistance  and 
inductance  matrices  of  the  patches,  thus  increeising  the  computation  time.  Since  we  use  a 
surface  integral  equation  method  for  high  frequencies,  we  do  not  change  the  distribution 
of  the  triangular  patches  for  the  cross-section  method  as  we  increase  the  frequency. 

For  high  frequencies,  we  use  a  coupled  surface  integral  equation  technique.  Under 
the  quasi-TEM  assumption,  the  frequency-dependent  resistance  and  inductance  result  from 
the  power  dissipation  and  magnetic  stored  energy,  which  can  be  czdculated  by  solving  a 
magnetoquasistatic  problem,  with  the  vector  potential  satisfying  Laplace’s  equation  in  the 
region  outside  all  the  conductors.  The  resistance  and  inductance  are  usually  given  by 
integrals  of  these  field  quantities  over  the  cross-sections  of  the  wires,  but  by  using  some 
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vector  identities  it  is  possible  to  convert  these  expressions  to  integrals  only  over  the  surfaces 
of  the  wires.  These  expressions  contain  only  the  current  at  the  surface  of  each  conductor, 
the  derivative  of  that  current  normal  to  the  surface,  smd  constants  of  the  vector  potential. 
A  coupled  integral  equation  is  then  derived  to  relate  these  quantities  through  Laplace’s 
equation  and  its  Green’s  function  outside  the  conductors  zmd  the  diffusion  equation  and  its 
Green’s  function  inside  the  conductors.  The  method  of  moments  with  pulse  basis  functions 
is  used  to  solve  the  integral  equations.  This  method  differs  from  previous  work  in  that  the 
calculation  of  resistance  and  inductance  is  based  on  power  dissipation  and  stored  magnetic 
energy,  rather  than  on  impedance  ratios.  It  will  therefore  be  more  easily  extended  to 
structures  where  non-TEM  propagation  can  occur. 

For  the  intermediate  frequency  range,  where  the  conductors  Me  on  the  order  of  the 
skin  depth,  were  found  it  very  efficient  to  interpolate  between  the  results  of  the  cross- 
section  and  surface  methods.  The  interpolation  function  was  based  on  the  average  size  of 
the  conductors,  measured  in  skin  depths,  and  was  of  the  form  1/(1  -f-  0.16a^/S^),  where  it 
a  is  the  average  cross-section  of  the  conductors,  and  S  is  the  skin  depth. 
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Coupled  Noise  Analysis  for  Ac^acent  Vias  in 
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Massachusetts  Institute  of  Technology 
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In  high-speed  digital  circuits,  hi^h  frequency  phenomena  affects  the  characteristics  of  the 
interconnections.  Physical  discontinuities  or  nonuniformities  in  the  connections  may  cause 
severe  reflections  when  they  can  no  longer  be  considered  as  conducting  wires,  but  behave 
as  transmission  lines  and/or  waveguides.  In  multilayered  digital  circuits,  vias  constitute 
one  of  the  most  commonly-used  class  of  interconnects. 

Vias  are  not  good  carriers  of  high-speed  signals.  They  cause  signal  distortion  and  re¬ 
flections  as  well  as  severe  degradation  in  the  high  frequency  components.  The  analysis  and 
modeling  of  a  single  via  and  some  quasi-static  or  quasi-TEM  analyses  of  single  via  config¬ 
urations  have  been  carried  out  previously.  To  date,  the  analysis  of  coupled  noise  between 
adjacent  vias  has  received  very  little  attention.  The  major  reasons  include  the  complexity 
of  multi-via  structures  and  the  difficulty  of  modeling  and  analyzing  them  accuratuy  over 
a  broad  frequency  range. 

The  electromagnetic  coupling  between  two  adjacent  vias  in  a  multilayered  integrated 
circuit  is  analyzed  by  means  of  equivalent  magnetic  frill  array  models  incorporated  with 
the  even-  and  odd-mode  approach.  Closed-form  expressions  for  the  coupled  noise  on  the 
passive  via  are  derived.  The  coupling  responses  in  the  frequency  domain  and  crosstalk 
waveforms  in  the  time  domain  for  some  multilayered  via  structures  are  calciUated  based 
on  these  formulas.  A  4-layer  experimental  model  is  constructed  and  measurements  are 
taken  for  the  transmission,  reflection,  and  coupling  responses.  The  measurements  show 
good  agreement  with  the  cidculated  results  over  a  frequency  range  of  up  to  18  GHz. 
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Electromametic  Radiation  from  Modules-on-Backplane 
Configurations  in  Conqiuter  Systems 


K.  Li,  M>  A.  Tu«ou4ji,  S.  Y.  Poh,  R.  T.  Shin,  and  J.  A.  Kong 
Research  Laboratory  of  Electronics 
Massachusetts  Institute  of  Technology 
Cambridge,  MA  03ISe 


Wiffh  Crequency  digital  signals  associated  with  computers  and  communications  equip¬ 
ment  nave  the  potential  for  electromagnetic  interference  (EMIl.  The  need  for  compliance 
with  the  radiated  emissions  limits  established  by  FCC  and  otner  regulatory  agencies  for 
such  equipment  have  led  to  increased  interest  in  understanding  and  minimising  EM  radia¬ 
tion.  As  dock  speeds  increase  beyond  100  MHs  and  edge-rates  fall  into  the  sub-nanosecond 
range,  the  EMC  enpneer  is  fac^  with  FCC  tests  for  frequendes  beyond  1  GHs  and  the 
prospect  of  increased  energy  for  frequendes  in  the  hundreds  of  megahertz  range. 

The  modules-on-backplane  configuration  is  a  coounon  assembly  scheme  for  computers. 
Embedded  in  these  modules  suid  bai^plane  are  reference  planes  which  could  be  as  large  as 
tens  of  centimeters.  The  reference  planes  configuration  is  of  interest  because  the  resonant 
frequendes  begin  in  the  low  hundreds  of  megahertz,  given  typical  dimensions.  The  voltage 
lev^  on  such  reference  planes  are  subject  to  noise  fluctuations  (reference  noise)  originating 
from  distiubances  such  as  inductive  discontinuities  at  connector  pins.  The  increase  in 
power  requirements  and  in  the  number  of  simultaneous  switching  drcuits  is  likely  to  result 
in  stronger  noise  sources.  A  simplified  model  is  used  to  analyze  the  modules-on-backplane 
coz^guration.  Conducting  planes  model  the  modules  and  the  backplane.  A  generator, 
which  is  positioned  at  the  module  and  backplane  interface,  models  the  inductive  voltage, 
as  a  noise  source,  developed  across  the  mod^e-backplane  connector. 

The  finitt-difference  time-domain  (FD-TD)  technique,  which  is  based  on  the  discretiza¬ 
tion  of  the  Maxwell’s  equations,  is  employed  in  this  analysis.  Nodal  electric  and  magnetic 
fidds  are  computed  on  a  rectangular  grid  within  a  computational  domain  subject  to  appro¬ 
priate  boundary  and  initial  conditions.  An  absorbing  boundary  condition  is  enforced  on 
the  outer  boundary  of  the  computational  domain  and  approximatdy  simulates  imbounded 
space.  In  the  implementation  of  the  FD-TD  technique,  the  areas  of  concern  addressed  in¬ 
dude  the  modeling  of  the  excitation  as  a  coaxial  line  probe  and  the  use  of  Prony’s  method 
to  obtain  late  time  responses  of  the  fidds  for  highly  resonant  configurations.  Radiation 
impedances  are  obtain^  from  the  voltages  and  currents  on  the  probe  and  radiation  pat¬ 
terns  are  calculated  using  Huygens’  prindple. 

Experimental  verification  is  msde  on  physical  structures  compatible  with  simulation 
models.  The  focus  is  on  corrdating  computed  resonant  frequendes  and  radiation  imped- 
^th  those  deduced  from  5-parameter  measurements  taken  on  a  network  analyzer. 
Hamation  impedances  as  functions  of  frequency  are  presented  for  various  configurations, 
wch  variations  indude  chan^ng  the  number  of  modules  and  the  dimensions  and  spadngs 
of  modules  and  backplane  and  the  addition  of  conducting  panels  to  modd  the  shidding  or 
ondosure  environment. 


study  provides  basic  information  pertaining  to  the  acceptability  of  noise  levds  in 
be  interconnect  environment  for  EMI-rdated  consideratioius.  Though  the  examination  of 
•  variety  of  configurations  it  will  also  affect  system-levd  pacloging  and  assembly  dedsions. 
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Analyses  of  Single  and  coupled  Dielectric  Bib  Wavegiudes 
and  Discontinuities  using  the  Finite-Difference  Time-Domain  Method 

Alex  Mou,  Y.  Eric  Yang,  Jin  A.  Kong 
Department  of  Electrical  Engineering  and  Computer  Science 
and  Research  Laboratory  of  Electronics 
Massachusetts  bistitute  of  Technology 
Cambridge,  MA  02139 

Millimeter  amd  submillimeter  transmission  structures  are  dominated  by  dielectric  waveg¬ 
uides.  At  these  frequency  bands,  conventional  microstrip  transmission  structures  such  as 
microstrip  line,  coplanar  line,  slot  line  and  finline  all  suffer  from  severe  ohmic  and  radia¬ 
tion  losses.  Analytic^  techniques  for  studying  dielectric  wavegmde  structures  include  the 
full-wave  integral  equation  formulations  (both  spectral  and  spatial  domain)  and  the  finite- 
difference  time-domain  (FDTD)  method.  Compared  to  the  integral  equation  methods,  the 
FDTD  method  is  more  straightforward  to  formulate  and  provides  an  alternative  path  to 
analyze  both  time  and  frequency  domain  response  in  the  presence  of  complex  structures 
such  as  discontinuities. 

In  this  paper,  electromagnetic  pulse  propagation  in  single  and  coupled  dielectric  rib 
waveguide  embedded  in  a  three-layer  stratified  ^electric  medium  is  simulated  and  analyzed 
using  the  FDTD  technique.  Time  domain  transient  response  of  a  modulated  gaussian 
pulse  is  observed  in  the  simulation.  We  first  examine  pulse  dispersion,  attenuation  and 
coupling  characteristics  for  the  transmission  line  structure.  Results  are  found  to  be  in  good 
agreement  with  those  obtained  from  integral  equation  methods.  We  also  extend  the  study 
to  an  optical  rib  waveguide  with  a  bend  discontinuity.  The  FDTD  simulation  is  used  to 
predict  the  electromagnetic  leakage  attributed  to  the  discontinuity. 

One  of  the  problems  facing  the  FDTD  technique  for  the  dielectric  waveguide  problem  is 
the  proper  absorbing  boundary  condition.  Existing  boundary  conditions  such  as  the  Mur, 
Enquist-Majda  condition  work  best  when  wave  incident  angle  is  normal  to  the  boimdary, 
which  is  generally  not  the  case  for  dielectric  waveguides.  A  comparison  between  different 
conditions  will  also  be  presented 
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The  use  of  equipment  enclosures  to  shield  against  emissions  and  external  energy  is  a 
common  practice.  Efforts  range  from  metal  impregnated  plastics  to  thick  steel  panels 
incorporating  waveguide  design  in  air  vents.  The  choice  of  complexity  depends  on  system 
performance  and  application.  The  question  has  often  been  raised  as  to  whether  high 
field  intensities  within  such  enclosures  will  affect  equipment  performance.  Clearly  this  is 
increasingly  pertinent  in  cases  of  heavy  shielding  where  resonator-like  structures  with  high 
quality  factors  ^Q)  result.  Moreover,  in  these  cases,  any  energy  lesJcage  may  be  highly 
frequency  selective  with  substantial  associated  field  strength. 

The  primary  motivation  for  studying  such  structures  is  that  computers  are  now  operat¬ 
ing  at  faster  speeds  and  consuming  more  power  resulting  in  significantly  increeised  levels  of 
power  at  higher  frequencies.  Also  typical  metallic  enclosure  configurations  have  dimeusions 
which  are  electrically  resonant  at  frequencies  in  the  hundreds  of  megahertz  range. 

This  paper  focuses  on  the  use  of  resistive  material  within  the  enclosure  walls  for  absorb¬ 
ing  some  of  the  contained  RF  energy,  thereby  avoiding  high  Q  enclosures  and  reducing  the 
likelihood  of  significant  field  buildup  within  the  enclosure.  The  energy  leakage  from  perfo¬ 
rations  on  the  walls  is  also  being  considered.  The  radiation  properties  of  these  structures 
are  observed  by  varying  the  number,  size  and  position  of  the  holes  in  the  enclosure. 

The  finite-difference  time-domain  (FD-TD)  technique  is  used  in  analyzing  this  problem 
and  the  metrics  are  total  radiated  power  and  field  patterns  within  the  enclosure.  The 
excitation  sources  are  line  current  source  for  the  two-dimensional  models  and  dipoles  for 
the  three-dimensional  cases.  A  Gaussian  amplitude  is  zissumed  in  order  to  obtain  multi¬ 
frequency  characteristics  and  to  exhibit  the  resonance  properties  of  the  model.  In  addition 
to  the  radiated  power  through  the  apertures,  the  power  absorbed  in  the  resistive  material 
is  also  calculated. 
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Abstract 


In  this  paper  a  new  method  bzised  on  the  wave  transmission  matrix  is  presented 
to  deal  with  nonuniform  multiconductor  line  system.  This  method  has  not  only  the 
ability  of  handling  nonuniformity,  but  also  the  advantage  of  handling  nonlinearity  effi¬ 
ciently.  The  wave  transmission  matrix  is  employed  to  characterize  the  multiconductor 
transmission  lines  and  the  telegrapher  equations,  therefore,  can  easily  be  turned  into 
algebraic  equations.  The  nonuniform  transmission  line  is  treated  as  many  small  seg¬ 
ments  of  uniform  transmission  lines  connected  in  cascade.  The  formulation  for  this 
method  is  described  in  detml.  Calculation  approaches  are  discussed.  Numerical  results 
are  presented  and  comparison  of  them  with  previous  work  confirms  the  validity  of  the 
method. 
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I.  Introduction 


In  very  large  scale  integrated  digital  circuits  and  systems,  multiconductor  transmis¬ 
sion  lines  are  commonly  used  as  signal  interconnections  between  chips,  chip  carriers, 
circuit  boards,  or  subsystems.  As  the  operating  speed  of  the  circuits  and  the  sys¬ 
tems  increases,  the  behavior  of  the  multiconductor  transmission  lines  employed  as  high 
speed  signaJ  carriers  will  have  significant  impact  on  the  performance  of  the  circuits  or 
the  systems.  A  multiconductor  transmission  line  usually  consists  of  cylindrical  and/or 
strip  conductors  embedded  in  an  inhomogeneous  medium  and  one  or  several  of  the 
conductors  serving  as  the  ground.  Inhomogeneous  transmission  lines  can  not  support 
pure  TEM  waves  but  quasi-TEM  waves  at  low  frequencies.  Even  if  the  medium  is 
homogeneous,  at  higher  frequencies  due  to  the  conducting  loss  and  higher  order  mode 
propagation,  the  lines  are  not  TEM  in  nature.  Since  a  rigorous  analysis  of  a  multicon¬ 
ductor  transmission  line  is  very  involved,  the  analyses  in  most  of  previous  works  and 
in  the  following  sections  are  based  on  the  quasi-TEM  approximation.  In  fact,  we  have 
seen  that  this  is  a  very  good  approximation  when  the  transverse  dimensions  of  the  line 
are  smtdl  as  compared  with  the  operating  wavelength  [1*11].  In  addition  to  simplifying 
the  analysis  of  the  multiconductor  lines,  the  quasi-TEM  approximation  is  also  desired 
for  handling  the  line  terminations  which  are  usually  described  by  circuit  quantities 
(voltages  and  currents)  rather  than  by  the  more  general  field  quantities  [12].  As  the 
frequency  of  operation  gets  higher,  the  modes  start  to  deviate  from  the  quasi-TEM 
modes  and  one  must  analyze  the  problem  with  rigorous  full-wave  approach  [13-15]. 

The  general  approach  on  the  analysis  of  the  multiconductor  line  is  as  follows:  uti¬ 
lizing  the  electromagnetic  field  theory  to  characterize  the  multiconductor  line  by  de¬ 
veloping  transmission  line  parameter  matrices,  such  as  the  capacitance,  inductance, 
resistance  and  conductance  matrices  per  unit  length  of  the  line,  and  to  establish  a 
mathematical  model,  i.e.,  a  set  of  equations  governing  the  behavior  of  the  line  system; 
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and  then  using  preferable  method,  for  example  the  modal  an3dysis,  to  solve  the  equa¬ 
tions  analytically  and/or  numerically  and  to  obtain  the  system  responses  in  the  time 
and  frequency  domains  through  corresponding  boundary  conditions.  Here,  we  shall 
emphasize  the  method  of  analyzing  the  responses  of  the  multiconductor  transmission 
line  in  the  frequency  domain  and  assume  that  the  parameter  matrices  of  the  line  are 
available  (refer  to  [16-21]. 

When  the  dominant  mode  in  the  multiconductor  transmission  lines  is  quasi-TEM, 
the  mathematical  model  which  fully  describes  the  line  systems  in  the  time  domain  or  in 
the  frequency  domain  is  the  well  known  telegrapher  equations  (or  refer  to  as  transmis¬ 
sion  line  equations).  For  uniform  multiconductor  transmission  lines,  the  corresponding 
telegrapher  equations  are  most  commonly  solved  by  using  the  so-called  modal  analysis 
in  the  time  domain  [1-7]  or  in  the  frequency  domain  [7-11].  However,  for  a  lossy  or 
frequency-dependent  line,  usually  the  modal  analysis  in  the  frequency  domain  will  be 
used  because  in  this  case  N  different  quasi-TEM  modes  propagating  on  a  transmission 
line  with  N  signal  conductors  and  one  ground  conductor  can  be  defined  in  the  fre¬ 
quency  domain  only.  In  the  modal  analyses,  the  kernel  of  the  approach  is  to  decouple 
the  telegrapher  equations  by  utilizing  the  method  of  characteristics,  and  then  based  on 
given  boundary  conditions  the  problem  can  be  solved  in  a  variety  of  ways,  for  example, 
by  solving  directly  [3,6-7],  or  by  network  analysis  methods  [10].  Other  approaches  for 
the  analysis  of  the  multiconductor  transmission  line  are  the  equivalent  circuit  model 
techniques  [22-23],  and  matrix  parameter  methods  (Green’s  function  methods)  [24- 
29].  The  modal  analyses  can  not  directly  deal  with  nonuniform  multi  conductor  line 
problems,  but  incorporating  with  other  techniques,  such  as  the  perturbation  method, 
iteration  method,  or  the  numerical  method,  it  is  still  possible  to  use  the  modal  analysis 
in  the  time  domain  for  some  nonuniform  systems  [3,5,30].  However,  the  more  efficient 
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methods  for  handling  the  nonuniform  transmission  lines  are  the  network  matrix  pa¬ 
rameter  methods  [24-29]  in  the  frequency  domain  and  the  spectral  method  [31]  in  the 
time  domain. 

We  shall  present  a  new  method  based  on  the  wave  transmission  matrix  to  deal  with 
the  nonuniform  multiconductor  line  system  in  this  paper.  This  method  has  not  only  the 
ability  of  hzmdling  nonuniformity,  but  also  the  advantage  of  handling  nonlinearity  in  the 
terminations  efficiently.  The  wave  transmission  matrix  is  employed  to  characterize  the 
multiconductor  transmission  lines  and  the  telegrapher  equations,  therefore,  can  easily 
be  turned  into  algebraic  equations.  The  nonuniform  transmission  line  is  treated  as 
many  small  segments  of  uniform  transmission  lines  connected  in  cziscade.  The  boundary 
conditions  are  formulated  with  linear  loads  in  spectral  domain  and  with  nonlinezir  loads 
in  time  domain. 

The  details  of  this  method  will  be  described  in  Sections  II.  The  treatment  of  the 
boundary  conditions  including  nonlinear  terminations  will  be  discussed  in  Section  III. 
CjJculation  approach  and  numerical  results  will  be  presented  in  Sections  IV  and  V. 


II.  Analysis  and  Formulation 


The  general  system  configuration  under  consideration  is  shown  in  figure  1,  The  sys¬ 
tem  consists  of  multiconductor  interconnection  lines  teiminated  with  lineax/nonlinear 
loads.  In  this  section,  we  use  the  wave  transmission  matrix  to  characterize  the  multicon¬ 
ductor  transmission  lines.  The  boundary  condition  or  termination  will  be  considered 
in  next  section. 

For  lossy  multiconductor  transmission  lines,  the  telegrapher  equations  in  the  fre¬ 
quency  domadn  can  be  expressed  as 


dz 


0  Z 
Y  0 


(1) 
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where  V  and  I  are  the  N  dimensional  voltage  and  current  vectors  on  the  N  transmission 
lines,  and  Z  and  Y  denote  the  N  x  N  impedance  and  admittance  matrices  per  unit 
length,  respectively.  The  matrices  Z  and  Y  are  symmetric.  For  uniform  multiconductor 
lines,  they  are  related  to  N  x  N  resistance  (R),  inductance  (L),  conductamce  (G),  and 
capacitance  (C)  real  matrices  per  unit  length  as  follows: 

Z  =  R  +  jwL  (2) 

Y  =  G  +  juC  (3) 

Now,  we  define  a  forward-going  and  a  backward-going  voltage  wave  vectors, 
and  V“,  as 

=  V  ±  ZqI  (4) 

or  the  voltage  V  and  the  current  I  as 

V=i(V++V-)  (5) 

I  =  -  V)  (6) 

where  Zq  is  defined  in  any  of  the  following  equivalent  forms: 

Zo  =  (ZY)^/2y-1  =  (ZY)"^/^Z  (7) 


Substituting  (5)  and  (6)  into  (1)  and  after  some  mathematical  manipulation,  we 
obtain  the  wave  trsinsmission  equations 

-5  r  v+ 1  ^  1  [  -(ZZq  1  +  ZqY)  (ZZo  ^  -  ZoY)  1  [  V+ 1  .8) 

dz[v-\  2  [-(ZZq^-ZoY)  (ZZo^-l-ZoY)]  LV 
It  is  easy  to  prove  that  the  diagonal  elements  of  the  2x2  partitioned  matrix  on  the 
right-hand  side  of  (8)  are  — (ZY)^/^  and  (ZY)^/^  ,  and  the  off  diagonal  elementary 
matrices  equal  zero.  Equation  (8)  can  be  rewritten  as 


d  [v+1  _  [-(ZY)^/^  0  fV'’’ 

dz  [V“J  “  [  0  (ZY)^/^]  [V“ 


(9) 
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where  0  denotes  a  matrix  with  all  the  elements  being  0. 

Prom  (9),  we  cam  see  that  the  forward-going  wave  vector  amd  the  backward¬ 
going  wave  vector  V“  a^e  decoupled  from  each  other.  This  is  only  true  for  uniform 
multiconductor  tramsmission  lines.  There  exists  coupling  between  these  two  wave  vec¬ 
tors  for  nonuniform  transmission  lines,  because  the  travelling  wave  in  either  direction  is 
continuously  reflected  by  the  nonuniformity  of  the  lines  along  the  propagation  direction. 

The  equations  in  (9)  are  first  order  derivative  matrix  equations,  it  is  obvious  that 
their  solutions  have  the  forms: 

V+(z)  =  V+.2:p(-K2)  (10) 

V-(z)  =  V^exp(Kz)  (11) 

where  Vq  are  constant  vectors  to  be  determined  by  boundary  conditions  of  the  lines, 
and  K  is  the  propagation  constant  matrix  of  the  line, 

K  =  (12) 

As  we  know  that  the  voltage  wave  vectors  at  any  two  positions,  zi  and  Z2,  on  the 
transmission  line  are  related  through  a  wave  transmission  matrix,  the  vectors  at  z, 
V^(z),  therefore,  can  be  obtained  from  the  wave  vectors  at  z  =  0,  V^(0),  in  terms  of 
the  following  expression 

■V+(r)' 

lv-(r)J 

where  [A,-^]  is  the  wave  transmission  matrix  of  the  transmission  line  segment  between 
z  and  0.  The  elementeiry  matrices  Ajy(i,j  =  1,2)  in  the  wave  treinsmission  matrix  can 
be  easily  obtained  by  utilizing  (10)  and  (11),  and  they  are,  respectively: 

A„(r)  =  czp(-Kz)  =  1  -  +  (K)^^ 

-(K)’^+  -  +  (-l)W^  +  -  (><“) 

Ai2(r)  =  A2i(z)  =  0  (14&) 


and 
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(14c) 


A22{z)  =  ezp(Kz)  =  1  + 

^  ^ - ^  ^ 

where  1  is  the  unit  matrix.  The  features  of  the  wave  transmission  matrix  are  very 
similar  to  the  chain  parameter  matrix,  those  identities  given  in  [24]  are  also  suitable 
for  the  wave  transmission  matrix. 

For  a  nonuniform  multiconductor  transmission  line,  we  can  approximate  it  as  many 
small  uniform  line  segments  with  different  characteristic  impedance  matrix  Zgjt  {k  = 
1,2, ...,m)  connected  in  cascade.  The  wave  transmission  matrix  [Ajj]  for  the  whole 
multiconductor  line  can  be  expressed  in  product  of  sub-wave-transmission-matrices 
[Ajt  jj]  characterizing  each  smadl  line  segment,  and  impedance  transition  submatrices 
ITjb.iyl, 


All  Ai2  _ 
A21  A22 

Am, 11  A,„^12  T-r  /rTjtii  Tjt,i2  Ajtu 
Am.21  Am  22  ,  U  Tjt  21  Tjt  22  Ajfe  21 

k=m-l  ’  •  j  1.  > 

Ak,n 

Ak,22_ 

)  (15) 

where  Ak,ij{i,j  = 

1,2)  have  the  same  forms  as  (14),  and  they  are: 

M.u  =  =  y; 

(16«) 

p=o 

Ak,\2  =  Ajfc^2i  =  0 

(166) 

and 

A,,,J  =  «p(K,A<,)  =  f;  i(K,A<,)f 

_ n  P' 

(16c) 

where  is  the  length  of  the  A:th  uniform  segment  at  z  =  zj^,  and  Kj^  is  the  propagation 

constant  matrix  of  the  fcth  segment  and  it  can  be  expressed  by  the  per  unit  length 
impedance  and  admittance  matrices  of  the  ibth  segment,  Zjt  and  Yjt,  as 

K,  =  (17) 
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In  (15),  the  elementary  matrices  Tjt^y  in  the  impedance  transition  submatrices  have 
the  following  forms: 

and 

Tfc.l2  =  21  =  2(^0.!:  -  Zo,l:+i)Zo_J  .  (19) 

where  Zo^jt  (^  =  1,2,...  ,m)  is  the  characteristic  impedance  matrix  of  the  l:th  uniform 
segment. 

Comparing  (15)  with  (13),  we  can  see  that  in  the  nonuniform  case  the  off  diagonal 
elementary  matrices  A12  and  A21  are  not  zero  matrices.  This  means  that  there  exists 
coupling  between  the  forward-going  and  backwaird-going  voltage  waves,  V'*’  and  V~ 
due  to  reflection  caused  by  the  nonuniformities  along  the  line.  The  coupling  between 
forward-  and  backward-going  waves  will  be  weak  if  the  nonuniformities  are  small. 

In  reality,  uniform  or  nonuniform  multiconductor  transmission  lines  are  usually 
terminated  with  linear  and/or  nonlinear  loads.  The  solution  for  the  multiconductor 
transnussion  line  systems  must  be  related  to  boundary  conditions  of  the  systems.  The 
approaches  for  treating  different  boundary  conditions  will  be  described  in  the  next 
section. 

III.  Boundary  Conditions 

In  the  previous  section  the  multiconductor  transmission  lines  were  characterized 
by  the  wave  transmission  matrix  A.  In  order  to  analyze  the  total  response  of  the 
system,  the  terminations  or  boundary  conditions  need  to  be  considered.  In  this  section 
we  formulate  the  boundaury  conditions  with  linear  loads  in  spectral  domain  and  with 
nonlinear  loads  in  time  domain. 
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We  assume  Thevenin  equivalent  voltage  sources  and  impedances  for  the  linear  ter¬ 
minations  of  the  multiconductor  transmission  line  system  of  length  I  as  shown  in  Figure 


1.  Thus  the  boundary  conditions  are: 

V(0)  =  Es  -  ZsI(O) 

(20) 

and 

V(/)  =  El  +  Zil(0 

(21) 

where  E5  and  denote  excitation  voltage  vectors  at  z  =  0  and  z  =  I,  respectively,  and 
Z5  and  axeNxN  diagonal  source  and  load  impedance  matrices  at  the  corresponding 
ends. 

For  linear  terminations,  we  can  determine  the  boundary  values  of  the  forward-  and 
backward-going  voltage  wave  vectors  from  (20)  and  (21)  using  (5),  (6)  and  (13).  The 
final  expressions  of  the  boundary  values  V'*'(0)  and  V“(£)  are: 

V+(0)  =  H  -  TsA^Ia,,  -  TsAjiil  +  rxAi2AJ,>)->riZo„ZJj‘AJj'i-‘ 

[TiEj  -  rsAJ,*(l  +  ri,A,2AJj‘)-'TiEi)l  (22) 

and 

v-(0  =  n  +  riAj2'A,2  -  riZo„zj,‘Aj2'(i  -  rsA2j‘A2,)-'r5Ajj'|-‘ 

(TiEi  -  riZo„Z  j,‘Aj2‘(l  -  rsAJ,*A2,)-‘TsE5)l  (23) 


where  and  are 

r*  7,  =  (1  +  Zs,iZ;-,>„  J->{1  -  Zs,iZJ>)  (24) 

and 

Ts,i  =  2(1  +  Zs,iZi,>o„)-’  (25) 

where  Zgj  and  Zom  the  characteristic  impedance  matrices  of  the  first  and  the  Izist 
line  segments.  When  deriving  (22)  and  (23),  the  following  identity  has  been  used. 

■A-I1A22  -  Ai2A22*A2iA22  =  ^Om^oi  (2®) 
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The  voltage  wave  vectors  at  any  point  z  on  the  line  can  be  obtained  from  (13)  and 


the  boundary  values  (22)  and  (23),  and  they  are: 


■v+(z) 

[v-(^)J 


Aii(^) 

A2i(^) 


Ai2 

A22 


V+(0) 

v-(o 


(27) 


where  matrix  [Ay(z)]  is 

fAiifz)  Ai2(«)]  _  I 

[A2i(z)  A22(«)J  ^"»i.22j  |mi=z/A/ 

n  ^*.12]  \^k,n  Afc  i2 

•••A  ^*,22]  [Aib,21  Ajfc_22 

ik=m|— 1 

(28) 


The  results  derived  here  are  in  frequency  domain  and  should  be  performed  for  each 
frequency  separately.  The  time  response  of  the  voltages  can  be  obtained  using  Fourier 
transform. 

For  nonlinear  terminations  of  the  miilticonductor  transmission  lines,  the  analysis  of 
boundary  conditions  utilizes  the  incident  wave,  and  reflected  wave,  which 

are  related  by  the  scattering  wave  matrix  S. 

=  SV*”*"  (29) 

or  in  terms  of  elementsiry  scattering  matrices  as, 

[V-/(0)]^rSH  Si2l[V‘-(0)]  (30) 

[v’-«/(/)J  IS21  S22J  [v‘"^(oJ 

where  ^*”*^(0)  and  V*”^(/)  are  the  incident  voltage  wave  vectors  at  z  =  0  and  z  =  /, 
and  (0)  and  (/)  are  the  reflected  voltage  wave  vectors  at  z  =  0  and  z  =  1. 
The  Sij  are  related  to  the  elementary  wave  transmission  matrices  Ay  as  follows 

511  =  -Aj2*A2i 

512  =  A22^ 

521  =  An  -  Ai2A22^A21 

522  =  Ai2A22^ 
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The  relation  between  the  reflected  and  incident  voltage  vector  in  time  domain  can 
be  written  as 

v*"*/  =  h  ♦  (31) 

where  h  is  the  inverse  Fourier  transform  of  the  scattering  matrix  S  and  asterisk  denotes 
convolution  with  respect  to  time. 

Consider  the  N  nonuniform  transmission  lines,  where  each  line  is  terminated  by  a 
voltage  source  and  a  load.  Using  the  relationship  of  the  incident  and  reflected  wave 
with  the  current  and  voltage,  the  following  expressions  can  be  obtained, 

v‘««  =  v'  +  c  -  z  •  i  (32) 

v"/  =  v'  +  e  +  z  *  i  (33) 


where  v'  is  the  voltage  vector  across  the  loads,  e  is  the  voltage  source  vector,  i  is  the 
current  vector  and  each  of  its  elements  is  a  function  of  the  corresponding  element  of  v'. 
The  matrix  z  is  the  inverse  Fourier  transform  of  the  diagonal  impedance  matrix  char¬ 
acterizing  the  terminal  impedance  lines.  The  first  N  diagonal  elements  of  z  correspond 
to  the  impedances  of  the  terminal  lines  connected  to  the  first  segments  of  the  trans¬ 
mission  lines  and  the  latter  N  elements  correspond  to  the  impedances  of  the  terminal 
lines  connected  to  the  last  segments. 

To  solve  the  equations  (31)-(33)  numerically,  we  discretize  the  time  variable  and 
convert  the  convolution  integrals  into  summation.  After  some  manipulation,  we  obtain 


ife-1  ife  ^ 

=  -*4  +  (1  -  hop*  {  -  [1  +  ho)  ^  I 


(34) 


m=0 


Tn=0 


In  the  above  equation,  the  subscripts  denote  the  discretized  time  step,  for  instant 
v'jt  =  v\kAt)y  where  At  is  the  time  step.  Equation  (34)  is  a  set  of  nonlinear  algebraic 
equations,  which  can  be  solved  by  means  of  numerical  techniques,  such  as  the  Newton 
Raphson  technique.  The  incident  voltage  used  in  equation  (34)  can  be  obtained  from 


k 

=  v'jt  +  Cjfe  -  ^ 
m=0 


(35) 
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The  above  formulation  would  also  work  for  linear  loads  for  which  the  current  will  be 
the  ratio  of  v'  to  the  load  impedance. 

rV.  Computational  Consideration 

The  major  part  in  the  numerical  implementation  of  this  method  consists  of  various 
matrix  operations  and  manipulations,  including  calculation  of  matrix  functions  such  as 
square-root  and  exponential  functions. 

In  the  calculation  of  wave  transmission  matrix,  the  square-root  manipulation  of 
a  matrix  is  involved  as  shown  in  (12).  A  formal  approach  to  carry  out  the  square- 
root  of  a  matrix  is  to  triinsform  the  matrix  into  a  diagonalized  matrix,  i.e.,  to  use  the 
characteristics  method.  Define  a  matrix  X  to  be  the  square-root  of  ZY,  then  we  have 

ZY  =  XX  (36) 

Assuming  that  the  eigenvzJue  and  eigenvector  of  ZY  are  and  D,-,  respectively,  they 
must  have  the  following  relationship 

D“^ZYD  =  diag(/??]  (37) 

where  D  is  a  matrix  consisting  of  the  eigenvectors  D,-,  i.e.,  D  =  [D], D2, . . . ,D^]. 
Substituting  (36)  into  (37),  we  obtain 

diag(^?]  =  D“^ZYD  =  D~^XDD"^XD  (38) 

From  above,  y/ZY  is  obtained  as 

X  =  =  Ddiag[/?,]D'^  (39) 

In  general,  a  matrix  function  can  be  defined  and  calculated  by  the  following  formula 

/(W)  =  Ddiag(/(A,),/(A2) . /(A;,)|D->  (40) 
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where  /(•)  denotes  a  function,  such  as  square  root  and/or  exponential  function,  and 
is  the  eigenvalue  and  D  is  the  matrix  consisting  of  eigenvectors  of  matrix  W. 

The  elementary  matrices  Ay  (*,j  =  1,2)  defined  by  (16a-c)  can  be  calculated 
directly  by  the  series  expansion.  The  matrix  series  are  absolutely  convergent  and 
their  convergence  rates  depend  upon  the  magnitude  of  eigenvalues  of  (*  = 

1,2,...,  N),  and  the  length  of  the  line  segments.  The  smaller  the  product  J 

the  higher  the  convergence  rate.  Therefore,  in  order  to  obtain  high  convergence  rate 
we  need  to  divide  the  uniform/nonuniform  line  into  lau^ge  number  of  short  segments, 
which  is  computationally  expensive.  The  exponential  matrix  functions  can  also  be 
czdcxilated  using  equation  (40),  as  in  the  calculation  of  square-root  function.  However 
the  latter  requires  the  calculation  of  the  eigenvectors  of  the  matrix,  a  procedure  very 
time  consuming  numerically. 

An  alternative  approach  to  calculate  the  elementary  matrices  is  to  make  use  of  the 
following  Sylvester’s  formula  (32).  U  &  N  x  N  matrix  W  has  N  distinct  eigenvalues,  a 
function  of  matrix  W  can  then  be  calculated  by 

N 

/(W)  =  fMZr  (41) 

r=l 

where  \r  are  eigenvalues  of  matrix  W,  and  Zr  are  N  x  N  matrices  given  by 

N 

Zr=U 

The  matrices  Zr  satisfy  the  following  relations: 

N 

Zr  =  1,  and  ZrZ,  =  0  (r  ,4  s)  (43) 

r=l 

The  approach  using  equation  (41)  avoids  the  calculation  of  the  eigenvector  and, 
since  it  uses  exact  formida  as  (40)  instead  of  series  expansion,  it  also  avoids  the  problem 
of  convergence.  Therefore  it  is  numerically  much  more  simpler  and  faster.  Moreover,  by 


W  -  A,1 

Ar-Ai 


(42) 
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using  (41)  in  the  caJculation  of  the  elementary  matrix  Ajj,  the  manipulation  of  square- 
root  and  exponential  of  matrix  can  be  combined  and  accomplished  in  one 

procedure.  For  the  cases  of  coupled  microstrip  lines,  the  eigenvalues  which  correspond 
to  coupled  modes  of  the  structures  are  always  non-degenerate  due  to  the  coupling  effect. 
In  our  numerical  implementation,  both  approaches  have  been  used,  if  the  eigenvalues 
of  the  matrix  product  ZY  are  distinct.  Same  results  are  obtained,  while  the- Sylvester’s 
formiila  signihcamtly  reduces  the  computation  time. 

To  validate  the  present  method  and  its  numerical  implementation,  here  we  consider 
a  pair  of  symmetric  nonuniformly  coupled  lines.  For  this  case,  the  elementary  matrices 
in  (16)  and  (17)  are  2x2  and  they  can  be  calculated  analyticadly.  Assuming  the  pair  of 
nonuniform  transmission  lines  have  the  following  impedance  matrix  Z  and  admittance 
matrix  Y : 


and 


Z  = 


Y  = 


'R{z) 

0 

t  4%  /  4 

Lm{z)' 

0 

R{z) 

T  yw 

Lmiz) 

L{z) 

fG(z) 

0 

C{z) 

Cmiz)] 

0 

G{z)\ 

-r  jw 

[Cm{z) 

C{z) 

(44a) 


(446) 


Then  the  square  root  of  ZY  can  be  obtained  by  the  following  expression 


K  =  V^=1 


7+  +  7- 
7-H  -7- 


7+  -7- 
7-1-  +  7- . 


where  7+  and  7-  are  the  eigenvalues  of  K  given  by 


(45) 


7+  =  +  ;■"(£(*)  +  £„(*))llC(z)  +  MC(.Z)  +  e„(r))]  (46a) 

7-  =  \/lH(r)  +  j<a(L(z)  -  L„(a)))lG(z)  +  iu.(C(x)  -  C„(z))|  (461) 

The  exponential  function  in  (16a,c)  are  pven  by  the  following  expression 

exp(-KAi)  =  -  f  exp(j+A£)  -  exp(^j-  Ai)  1  ^r) 

'  2  [exp(7-j-A^)  —  ezp(7_A£)  eip(7+A£) -f  eip(7- A£)  J  ^ 
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The  characteristic  impedance  matrix  Zq  will  be 


Zo  =  vWY"^ 

^  1  r 

2  -7-)  -  i^(7+  +7-)  2^^(7+  +7-)  -  ^(7+  -7-). 

(48) 

where 

AC  =  (G+ju,Cf+t^^Ci. 

For  this  example,  we  use  the  same  parameters  as  that  in  Fig.  8  in  [5].  The  total 
length  of  the  line  is  1cm.  The  elements  of  matrices  Z  and  Y  in  (44)  are: 

L(z)  =  5.73(1  —  ki(z))  (nH/cm) 

Lmiz)  =  ki{z)L{z)  (nH/cm) 

C{z)  =  1.81(1  +  M^))  (pF/cm) 

Cm{z)  =  -kc{z)C{z)  (pF/cm)  (49) 

kc{z)  =  0.2  —  0.05[1  —  cos(27rz)] 
ki^i)  =  0.2 
R{z)  =  G(z)  =  0 

A  pulse  source  with  amplitude  10  Volts,  both  rise  time  and  fall  time  of  10  ps,  and 
duration  20  ps  is  imposed  on  one  end  of  the  lines  and  all  the  ends  are  terminated  with 
50  fl  load.  The  time  waveforms  of  the  voltage  at  each  end  are  given  in  Fig.  2.  These 
results  compares  very  well  with  those  of  (5),  especially  for  the  near  end  and  fair  end 
coupling  noises.  The  numerical  calculations  of  the  matrix  are  compaired  with  analytical 
results  when  possible  and  very  good  agreement  is  found,  therefore  verified  the  numericaJ 
codes. 
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V.  Numerical  Examples 


In  this  section,  several  examples  of  nonuniform  multiconductor  microstrip  lines  with 
linear  or  nonlinear  terminations  are  considered.  In  the  analysis  of  any  transmission 
line  system,  we  need  first  to  construct  the  impedance  matrix  Z  and  admittance  matrix 
Y,  as  defined  in  (2)  and  (3)  by  resistance  (R),  inductance  (L),  conductance  (G), 
and  capacitance  (C)  matrices.  Many  methods  have  been  developed  to  calculate  the 
impedemce  and  admittance  of  microstrip  lines.  However,  the  formulation  developed 
in  the  previous  sections  does  not  depend  on  the  ways  those  parameter  matrices  are 
calculated.  Moreover,  by  introducing  different  effects  in  the  parameter  matrices,  the 
same  effects  can  be  incorporated  into  the  system  signal  responses.  For  example,  if  the 
inductance  and  capacitance  are  calculated  as  frequency  dependent,  then  the  dispersion 
can  be  included,  if  the  skin  effects  are  considered  in  the  czJculation  of  resistance  and 
conductance,  the  overall  losses  ais  a  function  of  frequency  can  be  determined.  In  the 
following  examples,  we  calculated  the  inductance  and  capacitance  in  a  way  similar  to 
that  of  [35].  The  resistance  and  conductance  are  calculated  by  the  strip  cross  section 
and  the  metal  conductivity.  The  elementary  matrices  Ay  axe  computed  by  utilizing 
(41)  (Sylvester’s  formula)  for  the  uniform/nonuniform  regions. 

As  first  example,  the  nonuniform  coupled  interconnection,  given  as  example  5  in 
[31],  is  analyzed.  A  ramp  function  incidence  with  a  rise  time  100  ps  and  an  amplitude 
1  volt  is  applied  on  the  active  line.  At  beginning,  we  employed  a  frequency  bandwidth 
from  DC  to  100  GHz,  and  obtained  the  time  waveforms  on  the  active  amd  passive 
lines  as  shown  in  Fig.  3(a)  and  3(b),  respectively.  Comparing  them  with  those  results 
presented  in  [31]  and  [33],  we  found  that  our  waveforms  are  sharper  than  theirs.  As  we 
use  less  frequency  bandwidth,  such  as  from  DC  to  15  Ghz,  the  resulting  waveforms  are 
smoother  as  shown  in  Fig.  4(a)  and  4(b),  and  they  aire  very  similar  to  those  given  in 
[31]  and  [33].  In  the  calculation,  we  use  80  small  uniform  segments  with  equal  length  to 
approximate  the  nonuniform  region  of  this  coupled  interconnection.  The  CPU  time  is  3 
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minutes  and  50  seconds  on  DEC  Station  3100  for  the  100  GHz  case  with  500  frequency 
points. 

As  discussed  in  Section  III,  multiple  line  system  with  nonlinear  terminations  can 
be  analyzed  efficiently  in  terms  of  scattering  matrix  and  the  corresponding  impulse 
responses  in  time  domain.  As  the  second  example,  we  consider  the  case  of  two  uniformly 
coupled  lines  with  nonlinear  terminations.  In  order  to  check  our  results,  we  used  the 
same  L,C,  R  and  G  parameter  matrices  and  nonlinear  load  as  that  given  in  [25].  Fig. 
5(a)  gives  the  near-end  and  far-end  voltage  responses  on  the  active  line,  and  Fig.  5(b) 
gives  those  on  the  parasitic  line.  Our  results  match  very  well  with  that  of  [25].  This 
further  confirms  that  our  method  is  efficient  and  valid  for  solving  multiconductor  line 
with  nonlinear  termination. 

In  the  following  example,  we  analyze  a  more  complicated  three  nonuniform  mi¬ 
crostrip  line  structure  with  both  lineu  and  nonlinear  terminations.  The  transmission 
line  system  and  the  cross  section  of  the  microstrip  are  depicted  in  Fig.  6.  The  thickness 
of  the  metal  strip  is  assumed  to  be  T  =  bfi/m.  The  surface  resistance  of  the  metal  strip 
is  R  =  30mCl/ square.  In  our  calculation,  the  1cm  long  nonuniform  microstrip  is  di¬ 
vided  into  50  segments  for  computing  the  elementwy  matrices.  The  incidence  voltage 
is  a  step  function  with  an  amplitude  of  1  volt.  First,  we  calculated  the  linear  load 
case  where  all  ends  are  terminated  with  50  fl  resistors.  Fig.  7(a)  shows  the  neair-ends 
response  voltages  and  Fig.  7(b)  shows  the  far-end  responses.  These  results  compare 
very  well  with  that  given  in  [34],  even  though  the  L,C,  R  and  G  matrices  are  calculated 
differently  and  the  strip  thickness  is  not  given  in  [34]. 

Next,  we  calculated  the  case  where  the  active  line  feu’-end  termination  was  replaced 
by  a  nonlinear  load  ,  but  the  rest  were  kept  the  S2ane.  The  nonlinear  load  is  a  diode 
with  nonlinear  characteristics  t  =  —  1),  where  /,  =  1mA  and  Vj’  =  25mV. 

In  Fig.  8(a)  and  8(b),  the  voltage  waveforms  on  the  fM  and  near  end  terminals  of  the 
three  lines  are  shown.  Comparing  Fig.  8  with  Fig.  7,  we  can  see  that  the  nonlinear 
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load  not  only  varies  the  waveform  at  port  4  on  the  active  line,  but  also  changes  the 
waveforms  at  the  other  ports. 

VI.  Conclusions 

A  new  method  based  on  the  wave  transmission  matrix  to  deal  with  nonuniform 
multiconductor  transmission  lines  has  been  developed.  One  of  the  advantages  of  this 
method  is  the  flexibility  of  handling  a  variety  of  nonuniform  multiline  structures,  be¬ 
cause  in  most  cases  it  is  possible  to  approximately  deal  with  a  nonuniform  multiline 
as  small  segments  of  uniform  multiline  connected  in  cascade,  and  the  total  wave  trans¬ 
mission  matrix  will  be  the  multiplication  of  all  the  submatrices  representing  each  small 
uniform  segment.  The  other  virtue  is  that  uniform  or  weak  nonuniform  multilines  will 
have  diagonalized  wave  transmission  matrices.  This  method  is  developed  on  the  vari¬ 
ables,  forward  going  and  backward  going  waves,  instead  of  the  voltage  and  current,  it 
is  consistent  with  the  method  solving  nonlinear  problems  in  terms  of  scattering  matrix 
and  the  corresponding  impulse  responses  (27,28].  They  can  be  easily  incorporated,  and 
then  an  approach  to  handle  nonuniform  multiconductor  transmission  line  with  nonlin¬ 
ear  termination  problems  can  be  achieved.  This  combination  method  is  very  efficient 
to  solve  the  nonuniform  structure  with  nonlinear  termination  problems,  even  if  the 
nonuniform  structure  is  a  lossless  or  very  low  lossy  system.  The  methods  based  on 
the  voltage  and  current  variables  will  have  heird  time  to  deal  with  the  very  low  lossy 
multiline  with  nonlinezu:  termination  problems. 

The  examples  given  in  the  previous  section  firmly  validate  that  this  method  is 
efficient  to  deal  with  nonuniform  multilines  and  is  significant  to  corporate  with  the 
method  solving  nonlinear  termination  problems.  It  is  possible  to  extend  the  combina¬ 
tion  method  for  solving  more  general  nonuniform  and  nonlinear  problems. 
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System  configuration  of  multiconductor  interconnection  lines  with  lin¬ 
ear/nonlinear  loads. 

Time  waveforms  of  transmitted  signal,  near  end  noise  and  far  end  noise 
on  the  coupled  lines  with  parameters  given  by  Equ.  (49). 

Time  waveforms  on  two  coupled  nonuniform  lines  with  frequency  band¬ 
width  from  DC  to  100  GHz,  (a)  active  line,  (b)  passive  line. 

Time  waveforms  on  two  coupled  nonuniform  lines  with  frequency  band¬ 
width  from  DC  to  15  GHz,  (a)  active  line,  (b)  passive  line. 

Time  waveforms  on  two  coupled  lines  with  nonlinear  terminations,  (a) 
active  line,  (b)  passive  line. 

Configuration  of  three  nonuniform  microstrip  system  with  linear  or  non¬ 
linear  loads. 

(a)  Near-end  and  (b)  far-end  response  voltages  with  linear  loads  for  the 
system  of  Figure  6. 

(a)  Near-end  and  (b)  far-end  response  voltages  with  nonlinear  loads  for 
the  system  of  Figure  6. 
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Abstract 

A  macroscopic  model  is  proposed  for  nonlinear  electromagnetic  phenomena  in  supriron- 
ductors.  Nonlinear  constitutive  relations  are  derived  by  modifying  the  linear  London's 
equations.  The  superelectron  number  density  as  a  function  of  applied  macroscopic  oirrent 
density,  ns(J),  is  derived  from  a  distribution  of  electron  velocities  at  a  certain  temperntmo 
T.  At  temperature  T  ^  OK,  the  function  rtg^J)  has  a  smooth  variation  near  the  macro¬ 
scopic  critical  current  density  Jc-  Agreement  has  been  found  between  this  ns(J,T)  mofh'l 
and  the  temperature  dependence  of  ris  in  the  two-fluid  model.  The  nonlinear  conductivities 
(Ts{J)  and  <Jn{J)  are  obtained  from  the  London’s  equation  with  the  modified  n.s(./)  func¬ 
tion.  Nonlinear  resistance  iZ(/),  kinetic  inductance  L^^I)  and  surface  impedance  Z*!/)  in 
thin  wire,  slab,  and  strip  geometries  are  calculated. 
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I.  Introduction 


Superconductors  have  great  potential  applications  in  many  fields.  For  example,  in  mi¬ 
crowave  integrated  circuits,  high-Q  resonators  and  microstrips  can  be  made  of  supercon¬ 
ductors  with  low  losses.  The  discovery  of  high-Tc  superconductors  has  also  made  a  liig 
impact  on  modeling  of  superconductors.  Problem  of  modeling  nonlinear  electromagnetic 
properties  of  a  superconductor  is  of  practical  importance  [1-9].  For  example,  in  application 
of  superconductors  to  high-Q  resonators,  fields  and  currents  are  very  large  at  resonance 
and  nonlinear  effects  are  inevitable. 

In  this  paper,  we  use  the  macroscopic  (classical)  theory  to  model  nonlinear  superconduc¬ 
tivity.  The  constitutive  relations  which  relate  electric  field  E  and  magnetic  field  D  to 
superconducting  current  density  Js  will  be  derived.  The  model  is  based  on  two  London 
equations  and  the  two-fluid  model.  New  problem  at  hand  is  to  incorporate  the  nonlinear 
effects  into  the  constitutive  relations.  Application  of  the  nonlinear  constitutive  relations 
will  provide  new  methods  for  studying  nonlinear  effects  in  superconductors. 

There  are  different  types  of  nonlinearity  in  superconductors.  One  type  of  nonlinearity  is 
displayed  in  polycrystalline  superconductors  where  granular  currents  are  involved  (1.2j. 
The  granular  nonlinearity  occurs  when  the  current  I  is  above  Ic  (here  Ic  is  the  threshold 
current  of  the  grain  junctions).  It  is  similar  to  a  p-n  junction’s  exponential  I  —  V  relation. 
It  is  also  found  in  type  II  superconductors  that  the  vortex  motion  in  the  mixing  state 
(between  the  superconducting  and  the  normal  states)  can  cause  a  nonlinear  V  ~  I  relation 
when  current  density  J  is  slightly  greater  than  the  critical  current  density  Jc  [3J.  The 
granular  or  vortex  nonlinearity  will  not  be  discussed  in  this  paper. 

Another  type  of  nonlinearity  is  intrinsic  for  edl  superconductors.  This  is  due  to  the  dc- 
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pendence  of  superelectron  number  density  n«  on  the  applied  current  density  J.  This 
nonlinearity  is  more  general.  There  have  been  few  papers  in  literature  directly  addressing 
this  nonlinearity  problem.  In  this  paper  we  will  focus  on  this  intrinsic  nonlinearity. 

It  is  necessary  to  point  out  that  this  paper  is  based  on  macroscopic  theory  with  certain 
hypotheses.  The  authors  do  not  intend  to  propose  a  rigorous  microscopic  theory  to  explain 
the  nonlinear  superconductivity.  The  purpose  of  this  paper  is  to  provide  engineers  with  a 
material-independent  macroscopic  model  for  nonlinearities  in  superconductors.  The  inoclel 
will  be  based  on  macroscopic  parameters.  The  model  is  intended  to  achieve  the  following; 
(1)  explain  the  nonlinear  voltage-current  (V  —  I)  relations  and  the  dependence  of  induc¬ 
tances  on  currents  (L  —  I  relations).  This  requires  derivation  of  a  complex  condtictivity 
d-(  J)  model;  (2)  explain  experimental  results  which  show  that  the  macroscopic  V  —  I  curve 
has  no  abrupt  transition  when  J  exceeds  J«.;  (3)  obtain  a  temperature  dependence  of  ng(T) 
which  agrees  with  what  is  assumed  in  the  two-fluid  model;  and  (4)  derive  a  general  scheme 
for  solving  nonlinear  electromagnetic  problems  in  superconductors. 

II.  Distribution  of  Electrons 

First  it  is  ^lssumed  macroscopically  that  when  an  electron  has  an  energy  E  <  Ef,  it  is 
in  superconducting  state  and  paired  with  another  electron;  when  E  >  Ec,  it  is  in  normal 
state.  This  characteristic  energy  Ec  is  a  function  of  temperature  T.  For  higher  T,  Er  is 
lower.  At  the  critical  temperature  T  =  Tc,  Ec  =  0.  Such  characteristic  Ec  resembles  the 
gap  parameter  A  in  the  BCS  microscopic  theory.  It  may  correspond  to  the  de-pairing 
energy  of  a  Cooper  pmr.  But  in  this  paper,  this  Ec  is  proposed  as  merely  an  assumption. 
The  origin  of  this  Ec  should  be  studied  in  microscopic  theory.  We  will  leave  this  problem 
and  use  Ec  as  a  starting  assumption. 
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We  assume  the  temperature  dependence  of  Ec  as 


£;c(T)  =  3.52fcflTc^l-~y,  {oxQ<T<Tc  (1) 

where  kg  is  the  Boltzmann  constant,  Tc  is  the  critical  temperature,  and  q  is  a  free  paranx- 
eter  which  may  depend  on  the  type  of  material  (a  >  0  is  required). 

In  m2icroscopic  theory,  the  energy  for  each  electron  consists  of  kinetic  energy  and  potential 
energy,  E  =  EfhEp  =  +  Ep,  where  m  is  the  mass  of  an  electron  and  v  is  its  veUaity. 

In  quantum  theory,  v  is  the  expectation  value  <tjf\v\ip>,  where  tl>  is  the  wavefunction  of 
the  electron  and  v  is  the  velocity  operator.  The  electron  energy  E  can  be  express<'(l  as  a 
function  of  temperature  T,  applied  current  density  J,  magnetic  field  if,  and  field  frequrnry 
f,E  =  E{T,  J,  /,  H).  Ec  is  the  origin  for  critical  vzJues  of  Tc,  Jc,  fc  and  He- 

Consider  a  case  where  the  electrons  have  one-dimensional  velocities,  for  example,  in  a  thin 
(radius  a  X  the  penetration  depth)  wire  within  which  the  current  flows  in  only  one 
direction.  The  electrons  in  the  sample  have  different  velocities  due  to  thermal  motion. 
The  average  velocity  of  all  electrons  is  nonzero  along  the  current  direction.  The  current 
density  J  is  a  macroscopic  quantity.  At  a  certain  point  r,  J  is  related  to  the  average 
velocity  of  the  electrons  in  a  small  volume  AV, 

j  =  ^  ,2) 

where  N  is  the  number  of  electrons  in  AF  and  q  is  the  charge  of  an  electron.  Not*-  that 
not  every  electron  hzis  the  same  velocity.  The  velocities  i;,-  of  the  electrons  obey  a  certain 
distribution.  Here  we  assume  that  the  number  of  electrons,  6N,  which  have  the  velocities 
between  v  and  Sv  obeys  a  one-dimensional  distribution.  If  we  define 


(31 


then  (10,11] 


n(u)  =  n(t;,t;yi,no,m,fcflr)  (4) 

where  Vyi  is  the  average  velocity  of  the  electrons,  m  is  the  mass  of  an  electron,  and  no  is  the 
total  number  density  of  electrons  in  AK  Figure  1  shows  a  possible  n(v)  function  where  n{v) 
is  assumed  to  be  a  continuous  distribution  of  a  Gaussian  form.  The  specific  mathematical 
form  of  the  distribution  n(w)  may  be  obtained  from  the  microscopic  quantum  theory.  If 
the  wavefunction  of  each  electron  0,-  is  taken  into  account  for  a  distribution  n(e,  ),  where 
Vi  —<  0i|t;j0j  >,  such  formula  can  be  derived  from  statistical  physics.  The  numerical 
results  shown  in  this  paper  are  obtained  by  using  the  distribution  in  Figure  1. 

We  can  relate  the  two  macroscopic  quantities,  J  and  by 

/+00 

qn(v)vdv  =  gnQUyi  (5) 

-oo 

at  AV. 

The  T  dependence  of  n(u)  requires  that  at  T  =  0, 

n{v)  =  nQ6{v),  (G) 

a  delta  function  (12).  For  higher  T,  the  spread  of  the  distribution  is  wider. 

III.  ns{J)  and  Critical  Current  Density  Jc{T) 

Since  the  electrons  have  different  velocities,  and  therefore  different  energies  at  a  certain 
applied  current  density  J,  the  electrons  will  not  all  exceed  the  Ec  at  the  same  time.  The 
number  density  of  superelectrw..s,  ns,  does  not  disappear  abruptly  when  J  exceeds  th  ' 
macroscopic  critical  current  density  Jc-  Only  when  T  =  0  and  n(v)  becomes  no6(t’).  ^si-H 
shows  an  abrupt  drop  to  zero  at  a  critical  JdT  =  0):  for  J  <  7c(0),  ns  =  ng;  and  for 
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J  >  Jc(0),  ris  =  0.  In  experiments,  since  absolute  T  =  OK  is  not  achievable,  a  completely 
sharp  transition  of  n*  is  not  observed. 


The  JciT)  is  related  to  the  average  velocity  and  characteristic  energy  Ec.  First,  we 
define  a  characteristic  velocity  Vc  for  a  single  electron  corresponding  to  Ec-  If  the  potential 
energy  is  included  in  the  case  of  maximum  kinetic  energy,  ^tnv^  =  Ec-  Hence,  jv  = 
y/2Ec/rn-  Since  J  =  nQqvy^,  we  define 


Jc{T)  =  qriQVc  =  qriQ 


J 


IMkgTc 

m 


(7) 


as  the  critical  current  density.  This  Jc  is  different  from  the  Ginzburg- Landau’s  de-pairing 
current  density.  Here  JdT)  is  a  quantity  derived  from  the  hypothetical  characteristic 
energy  Ec  via  the  classical  kinetic  energy  expression.  For  T  ^  0,  ns{J)  is  a  smoothly 
varying  function  and  the  sample  is  partially  superconducting.  The  Jc(T)  has  a  different 
characteristics  from  J<r  at  T  =  0.  This  smooth  varying  feature  of  the  /  -  K  curve  has  been 
widely  observed  in  experiments  (5.G]. 


Since  VjH  =  J/qnQ,  the  number  density  of  superelectrons  (here  we  count  the  single  electron 
density  rather  than  the  pair  density)  ns{J)  can  be  derived  from  the  velocity  distribution. 
The  superelectrons  are  those  whose  energy  are  lower  than  Ec-  In  the  one-dimensional  ca.se. 

rAT) 

ndJ)=  I  n{v)dv  (S) 

J-vAT) 

In  general,  this  integral  can  only  be  evaluated  numerically. 


At  T  =  Tc,  Vc  =  0,  hence  =  0  for  all  J's.  By  using  the  electron  velocity  distribution  in 
Figure  1,  ns(J)  is  plotted  in  Figure  2  for  four  different  temperatures.  Here  J  is  normalized 
by  the  Jc  at  T  =  0.  At  T  =  0,  n.,  =  ng  for  J  <  Jc(0).  This  is  what  has  been  predicted  I)y 
the  two- fluid  model. 
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The  temperature  dependence  of  n*  assumed  in  the  two-fluid  model  is  that 

n.m  =  no  [l -(!•)<  . 

■‘C  . 

The  ris(T)  function  is  compared  with  our  distribution  model  in  Figure  3.  For  a  fixed  .7, 
the  rtg(T)  curve  matches  the  two-fluid  model  for  a  =  3/2. 


IV.  Conductivities  <Ts(J)  and  an(J) 

To  derive  the  constitutive  relations  of  superconductors  for  electromagnetic  fields,  we  will 
use  the  London’s  equations.  The  London’s  equations  are  derived  from  the  fundamental 
Newtonian  dynamics  and  the  Meissner  effect.  They  can  also  be  derived  from  quantum 
mechanics  by  introducing  a  canonical  momentum  [3,8].  If  we  do  not  consider  the  Lorentz 
force  due  to  magnetic  field,  the  linear  London’s  equations  are  valid  from  the  Drude  mod<'l 
and  the  Newton’s  second  law. 


V  X  Js  =  (10) 

A 

where  E  and  D  are  the  total  electric  and  magnetic  fields,  respectively,  Jg  is  the  cunrnt 
density  due  to  superelcctrons,  and  A  =  rng/q^ris-  The  subscript  s  denotes  that  the  quantity 
is  of  superelectrons.  The  nonlinearity  will  be  included  in  A(/)  = 


A(/)  = 


Once  ns(J)  is  known,  A(7)  can  be  derived.  Hence,  the  nonlinear  constitutive  relations  arc 
obtained.  Nonlinear  effects  come  in  A  and  n,. 


Substituting  eq.  (8)  of  risiJ)  in  A,  we  derive  the  A(  J)  function,  which  is  plotted  in  Figure 
4.  The  penetration  depth  goes  to  infinity  when  J  S>  JdT).  Conductivities  crs(J)  aiKl 
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(12) 


(Tn{J)  are  derived  from  the  London  equations  for  time-harmonic  (e  fields 

iq^n, 

(jg  — 

urrig 


and 


On  = 


q^rinT 
m(l  —  iu)T) 


(13) 


where 


n{y)dv 


and  the  number  density  of  normal  electrons  nn{J)  =  >^o~”s(*^)  from  conservation  of  charge. 
The  above  equations  are  the  main  results  of  this  paper.  The  conductivities  are  plotted  in 
Figures  5a  and  5b,  where  r  is  the  transport  time  or  the  mean  scattering  time  of  normal 
electrons.  (Ts(J)  and  CniJ)  reflect  the  current  dependence  of  Ug  and  nn,  respectively.  These 
plots  are  for  T  =  88  K  very  close  to  the  critical  temperature  Tc  =  90  K.  The  nonlinear 
behavior  is  easy  to  see.  Total  conductivity  is  a  complex  number 


d(J)  =  a„iJ)  -f-  cTsiJ) 


(14) 


V.  Geometry  Effect  in  Nonlinear  Relations 

1.  Thin  wire 

For  a  thin  wire  superconductor  with  radius  a  A,  the  current  density  J  can  be  assumed 
independent  of  the  radius  p.  The  thin  wire  is  an  ideal  case  since  a  A  is  not  practical. 
Usually  r  w  10“^^  sec,  and  u  <  10^  Hz,  hence  wt  <  1  which  is  the  quasi-static  case. 
an  w  ^  .  Here  we  can  not  assume  ^  |<7n(  since  we  are  studying  the  transition 
region  where  rig  may  become  very  small. 


8 


(15) 


Resistance  per  unit  length  of  the  wire  can  be  written  as 

o  _  _ l^nl 

IknP  +  \(Tsfi]ira^' 

The  kinetic  inductance  per  unit  length  of  the  wire  can  be  written  as 

r  k<l 

^  (kn|2  +  k*l^]u?7ra2‘ 

In  this  case,  the  intemsd  and  external  inductances,  and  Lex,  which  ^lre  related  to  the 
energy  stored  in  the  magnetic  field,  are  independent  of  I  and  therefore  are  linear.  Figure 
6  shows  the  calculated  R  —  I  and  Lj^  —  I  curves  from  eqs.  (15)  and  (16),  where  ut  =  10~^, 
I  =  Jna?.  The  R  —  J  curve  appears  very  nonlinear  because  T  =  88  K  is  very  close  to 
Tc  =  90  K.  One  interesting  behavior  in  the  Lj.  —  J  curve  is  that  the  maximum  Ls  appears 
near  Jc-  This  can  be  explained  from  eq.  (16).  At  J  Jc(T),  ris  ^  rinuir,  hence  Lf.  l/ug 
is  small.  At  J  ^  Jc(T),  nnujr  ';$>  rig,  hence,  L;.  «  Tig/(nn)‘^  is  also  small.  When  rinLjT  and 
Tig  are  comparable,  a  maximum  of  may  be  achieved. 


A  special  case  is  at  DC  when  u;  =  0.  FVom  eq.  (12),  erg  — »  oo.  Trom  the  first  London 
equation,  E  =  'gj-  =  0.  When  J  ^  Jc,  rtg  —*  0,  then  it  is  possible  for  E  0.  If  has 
a  small  time  fluctuation,  then  ^  is  not  exactly  zero  and  the  above  quasi-static  disr\ission 
can  be  applied. 


2.  Slab  and  thin  film  strip 


The  above  discussion  has  not  taken  into  account  the  geometry  of  the  superconductor.  For 
superconductors  with  finite  dimensions,  e.g.,  a  slab  in  Figure  7,  the  non-uniform  distribu¬ 
tion  of  current  density  Jy{x)  will  cause  the  nonlinearity  to  occur  earlier  in  the  R  —  I  and 
Lk  —  I  curves.  The  internal  and  external  inductances  and  Lex  will  also  be  nonlinear 
since  Jy{x)  is  determined  by  I, 


~  II  ~  I  J  —  J  j  dxa{x)E{x)  (17) 
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where  B(x)  is  the  electric  field.  For  a  uniform  current  density,  Jy(x)  =  7^  =  I/Azd 

where  d  is  the  thickness  of  the  slab.  For  non-uniform  Jy{x),  Jy  is  bigger  thain  at  the 
edge  of  the  slab  and  ns{J)  is  smaller  at  the  edge.  Therefore,  for  the  same  magnitude  of  /, 
inhomogeneous  Jy{x)  will  exceed  Jc  at  some  I’s  even  when  Jj\  is  less  than  Jc- 

If  a  three-dimensional  velocity  distribution  is  considered  for  a  current  flow  in  the  y  direc¬ 
tion,  Jy  ~  and 

^im^nCiJ)  =  nQ6{vx)6{vy  -  (IS) 


ns(Jy)  =  I  dvy  I  _ dvx  I  _ n{v)dvz  (19) 

J-Vc  J-y/vl-vl  J -  s/v‘{-vl-v\ 

with  =  Jy/qriQ. 

Under  the  quasi-static  condition  wr  <C  1,  if  a  current  I  is  applied  along  y  direction  in 

the  slab,  a  magnetic  field  Hz{x)  will  result  from  the  applied  current  Jy{x).  Note  that  in 

the  finite-width  slab  of  a  linear  superconductor,  even  at  u;  =  0,  current  Jy{x)  has  a  non- 
uniform  distribution,  due  to  the  penetration  depth  A.  From  the  second  London  equation 
(10)  and  VXi4  =  B,  Vx//  =  J,  under  the  London  gauge,  V  •  A  =  0  and 

A  =  -\Js,  (20) 

V^[7sA2(7)|  =  7, =  7  (21) 

where  A  is  the  vector  potential.  In  general,  this  equation  can  not  be  solved  in  a  closed-form. 


We  first  consider  the  solution  of  a  linear  superconductor  for  n®  rinUT. 


and 


Jy{x] 


I  cosh(r/A) 
2Asinh(d/2A) 


for  |x|  <  (f/2 


Hz{x) 


I  sinh(z/A) 
2sinh(7/2A) 


for  |i|  <  7/2. 


(22) 


(23) 


An  iterative  scheme  for  solving  the  nonhnear  problem  is  to  initially  use  the  linear  solution 
Jy{x)  in  eq.  (22)  to  calculate  ns{J)  in  eq.  (8)  and  A(x)  in  eq.  (11).  Then  A  in  eq.  (22)  is 
substituted  with  the  A(x)  obtained  to  calculate  Jy{x).  Then  eq.  (8)  is  used  to  calculate 
ns(Jy(i)]  again.  This  procedure  is  repeated  until  Jy{x)  converges.  Figure  8  shows  the 
results  of  J(x)  and  <r4(x)  obtained  from  the  iterative  method. 

The  impedance  of  the  slab  is  defined  as 

z  =  -  (24) 

and  resistance  R  =  Re{Z)  and  kinetic  inductance  L|.  =  —lTn{Z)/u.  Figure  9  shows  the 
R  —  I  and  —  I  relations  from  the  final  convergent  Jy{x). 


For  different  thickness  d,  the  nonlinear  curve  will  be  different.  We  have  found  that  for 
smaller  d/A,  R  —  I  and  —  I  relations  are  more  norxlinear.  This  is  understandable  since 
Jy  is  bigger  for  smaller  d  for  a  given  applied  current  I. 


Surface  impedance  of  a  superconductor  is  defined  as 


where 


c  =  c  + 


id{J) 


Surface  resistance  is  Rs  =  Rc{Zs)j  and  surface  inductance  is  Ls  =  —Im{Zs)l^.  For  wr  = 
10~^,  Rs  and  Ls  are  plotted  in  Figure  10.  The  nonlinear  region  appezu-s  near  Jc  where  the 
sample  is  partially  superconducting.  In  the  normal  or  superconducting  states,  the  relations 
are  linear.  The  Ls  is  different  at  normal  and  superconducting  regions.  Experimental 
results  [13]  are  compared  in  Figure  11.  As  it  shows,  this  model  agrees  .with  the  trend 
of  the  experimental  results.  Since  some  parameters  (o,  r,  etc.)  are  material  dependent. 
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by  adjusting  these  parameters,  a  better  fit  between  the  theory  and  measurement  may  be 
found.  Surface  impedance  Zs  is  a  good  description  of  superconductors  since  most  of  the 
currents  and  fields  are  confined  within  the  penetration  depth  A  from  the  surface. 


Another  case  of  interest  is  at  very  high  frequencies  and  On  for  T  >  Tc  is  not  very  big 
(which  is  true  for  the  ceramic  high  Tc  superconductors)  so  that  e  can  not  be  neglected  in 
eq.  (26)  for  Jy  near  Jc-  Since  ^  ^  0  in  the  Maxwell’s  equations,  a  wave  equation  has  to 
be  considered.  We  will  solve  the  guided  wave  case  where  the  electric  field  Ey  is  decaying 
away  outside  the  slab. 

The  wave  equation  for  E  is 

+  jtgcr-B  =  0  (27) 


where 


Cr 


=  1  + 


icr{J) 

weo 


(2S) 


If  we  replace  E  by  aJy  in  eq.  (27),  (27)  will  look  very  similar  to  (21).  We  can  also  u.se 
an  iterative  scheme  to  solve  for  the  nonlinear  a  case.  First,  we  assume  cr  is  indeix'iideiit 
of  J  and  solve  for  Jy{x).  Second,  we  calculate  <7(z)  =  (7[7y(i)].  Third,  we  solve  tlu-  wave 
equation  for  the  inhomogeneous  medium  problem  and  obtain  Jy(i)-  These  steps  will  lx- 
repeated  until  Jy{x)  and  ^•(.t)  converge.  This  is  a  feed-back  process.  At  the  edge,  increase 
in  Jy  will  cause  n*  to  decrease,  which  causes  a®,  and  hence  Jy  to  decrease.  This  will 
continue  until  the  a  stable  solution  is  obtained. 


This  procedure  is  similar  to  solving  coupled  two  differential  equations  in  the  GL  theory 
where  the  equation  for  A  is  essentially  the  same  £is  the  wave  equation  (27)  for  E  or  Jy. 
The  difference  is  the  second  equation  for  where  rp{A)  (ns(  J)  )  is  to  be  derived.  Wc-  have 
derived  ns{J)  in  this  paper  from  the  velocity  distribution  assumption. 
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The  results  for  the  slab  geometry  can  be  used  to  study  a  microstrip  geometry.  If  tlie 
dimension  of  z  in  Figure  7  is  reduced  to  <f  A,  the  current  will  still  be  uniform  in  the  r 
direction.  Therefore,  the  same  results  can  be  applied. 

This  nonlinear  model  can  also  be  applied  to  study  the  dependence  of  R  and  L  on  the 
magnetic  field  H.  Once  the  relation  between  the  magnetic  field  H  and  the  current  J  is 
determined,  He  and  Jc  can  be  related  and  the  above  discussion  is  directly  applicable. 

Although  the  model  presented  is  classical,  the  corresponding  quantum  statistical  distrilju- 
tion  can  be  used  to  derive  the  velocity  (energy)  distribution.  Discrete  distribution  may  be 
needed  if  the  energy  is  quantized. 

VI.  Conclusion 

A  macroscopic  model  is  proposed  for  nonlinear  constitutive  relations  in  supercondm  tors. 
Distribution  of  electron  velocities  is  used  to  derive  the  dependence  of  superelectron  df'nsity 
Tig  on  applied  macroscopic  current  density  J.  Complex  a(  J)  is  obtained.  The  geometry  of  a 
superconductor  will  introdtice  non-uniform  distribution  of  the  current  density  J.  Therefore, 
nonlinearity  will  be  enhanced  at  the  edges  and  surfaces  of  superconductors.  Using  this 
macroscopic  model,  a  solution  scheme  for  electromagnetic  properties  of  superconductors 
has  been  proposed. 
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Figure  Captions 


Fig.  1.  An  assumed  distribution  of  electron  velocities  for  T  =  20  K. 

Fig.  2.  Superelectron  number  density  n,(  J)/no  at  different  temperatures.  In  this  example, 
a  =  3/2,  zuid  Tc  =  90K.  Jc(0)  is  from  eq.  (7). 

Fig.  3.  Comparison  between  the  two-fluid  model  and  this  model  for  ns(T)/rio  at  .7  =  0. 
Here  a  =  3/2  and  Tc  =  90  K. 

Fig.  4.  Penetration  depth  A  as  a  function  of  current  density  J  T  =  88K.  Here  Tc  =  90 
K. 

Fig.  5.  At  T  =  88K,  Tc  =  90K,  (a)  super-conductivity  ag  as  a  function  of  J.  (b)  nonnal 
conductivity  <rn  as  a  function  of  J. 

Fig.  6.  At  T  =  88K,  Tc  =  90K,  (a)  resistance  of  a  thin  wire  R  as  a  function  of  .7.  (I)) 
kinetic  inductance  of  a  thin  wire  as  a  function  of  J. 

Fig.  7.  A  superconducting  slab  with  a  thickness  d. 

Fig.  8.  At  T  =  80K,  Tc  =  90K,  d/X{J  =  0)  =  2,  (a)  current  density  Jy{x)  distribution  in 
a  slab,  (b)  super-conductivity  <rs(i)  distribution  in  a  slab. 

Fig.  9.  For  d/\(J  =  0)  =  0.5, 1,2,  7^  =  90K,  (a)  resistance  of  a  slab,  72,  as  a  function  of 
applied  current  intensity  /  at  T  =  80K;  (b)  kinetic  inductance  Lj^  of  a  slab  as  a  binction 
of  applied  current  intensity  /  at  T  =  80K.  Here  Tc  =  90K. 

Fig.  10.  At  T  =  8SK,  Tc  =  90K,  (a)  Surface  resistance  72s  as  a  function  of  J,  (b)  surface 


inductance  £«  as  a  function  of  J . 


Fig.  11  as  a  function  of  surface  current  density  A"  at  T  =  77K.  The  circles  are  the 
measured  data  from  reference  [13].  The  sample  is  a  YBCO  at  /  =  10.4GHz  with  Tc  =  92K. 
The  solid  curve  is  from  this  model,  where  a  =  3/2,  and  ut  =  2.8  x  10“^.  The  surface 
current  density  is  calculated  from  J  with  a  thickness  of  A(7). 
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